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A Numerical Simulator for Three-Dimensional 
Unsteady Flows through Vibrating Blade Rows 


Summary 


The three-dimensional, multi-stage, unsteady, turbomachinery analysis, TURBO, has 
been extended to predict the aeroelastic and aeroacoustic response behaviors of a single 
blade row operating within a cylindrical annular duct. In particular, a blade vibration 
capability has been incorporated so that the TURBO analysis can be applied over a solution 
domain that deforms with a vibratory blade motion. Also, unsteady far-field conditions have 
been implemented to render the computational boundaries at inlet and exit transparent to 
outgoing uns teady disturbances. The modified TURBO analysis is applied herein to predict 
unsteady subsonic and transonic flows. The intent is to partially validate this nonlinear 
analysis for blade flutter applications via numerical results for benchmark unsteady flows, 
and to demonstrate the analysis for a realistic fan rotor. For these purposes, we have 
considered unsteady subsonic flows through a 3D version of the 10th Standard Cascade, 
and unsteady transonic flows through the first stage rotor of the NASA Lewis, Rotor 67, 
two-stage fan. 
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1. Introduction 


The development of analyses to predict unsteady flows through turbomachinery blade 
rows has been motivated primarily by the need to predict the aeroelastic (flutter and forced 
vibration) and aeroacoustic (sound generation and propagation) characteristics of the blad- 
ing. Accurate and efficient aerodynamic analyses are needed to determine the unsteady 
loads that act on the blades and the unsteady pressure responses that persist upstream and 
downstream of the blade row, for various sources of unsteady excitation. The latter include 
structural (blade) motions and aerodynamic disturbances at inlet and exit that carry energy 
towards the blade row. 

The computational resources required to simulate nonlinear unsteady flows continue to 
prohibit the use of such simulations in detailed aeroelastic or aeroacoustic design studies. 
Thus, for the most part, the unsteady aerodynamic analyses that are being used in tur- 
bomachinery aeroelastic and aeroacoustic design prediction systems are based on linearized 
inviscid flow theory [Ver93], which has evolved to the point that three-dimensional, lin- 
earized, Euler analyses are being developed [HL93, HCL94, Sre96, MV97]. Linear analyses 
meet the needs of turbomachinery designs for efficient unsteady aerodynamic response pre- 
dictions. However, of necessity, such analyses ignore potentially important physical features 
of unsteady flows, including the effects of moderate to large amplitude unsteady excitation 
and the effects of viscous-layer displacement and separation. 

Time-accurate, nonlinear, Euler and Reynolds-averaged Navier-Stokes analyses, are there- 
fore needed to understand and predict the relative importance of nonlinear and viscous ef- 
fects on the unsteady flows associated with blade vibration and blade-row noise generation. 
Since the mid 1980’s, a number of such analyses have been developed for turbomachin- 
ery configurations. These have been applied to predict flows through single blade rows in 
which the unsteadiness is caused by prescribed blade vibrations [HR89, Sid91, HD93, GV94, 
PGW96, GC96, AV96, BSK97] or by prescribed aerodynamic disturbances at the inflow 
or outflow boundaries [Gil88, DV94, CCA94], and flows through aerodynamically coupled 
arrays in which the unsteadiness is caused by the relative motions of adjacent blade rows 
[Rai87, Rai89, JW89, JW90, JHW92, CW93]. 

These recent and important advances in the numerical simulation of unsteady flows 
demonstrate the power and potential usefulness of nonlinear unsteady aerodynamic codes. 
Most of the related activity to date has been focused on developing numerical strategies for 
solving the Euler and Navier-Stokes equations and on implementing these strategies into use- 
ful codes. Although there is still a need for improvements in algorithm speed and accuracy 
and, more importantly, in the treatment of flow boundary conditions, a major focus of con- 
tinuing work must be placed on validating the capabilities of modem Euler/Navier-Stokes 
solution schemes for accurately predicting turbomachinery flow phenomena. Once validated, 
Euler and Navier-Stokes analyses for turbomachinery unsteady flows can provide engineers 
with useful insights into the impact of nonlinear and viscous effects on blade vibration and 
discrete-tone noise generation. These analyses would also provide a test-bed for evaluating 
and improving the linearized models that are being developed for use in aeroelastic and 
aeroacoustic design prediction systems. 

Under the present effort, we have modified and applied the nonlinear unsteady analysis, 
TURBO, to predict unsteady flows through single vibrating blade rows. TURBO is a multi- 


3 



stage turbomachinery code that has been constructed, as part of a long range research effort 
[Jan89, JW89, JW90, JHW92, CW93, CCA94] conducted at Mississippi State University 
(MSU), to simulate the complex unsteady flow phenomena occurring in turbomachines. 
The Euler, the thin-layer Reynolds-averaged Navier-Stokes equations or the full Reynolds- 
averaged Navier-Stokes equations axe solved using an implicit, cell-centered finite volume 
scheme, in which inviscid or convective flux Jacobians are evaluated using flux vector splitting 
and inviscid residual fluxes are evaluated using Roe’s [Roe81] flux difference splitting to 
form a higher-order TVD scheme. Viscous or diffusive fluxes can be treated either explicitly 
[CW93] or implicitly [CCA94]. Newton subiterations are used as part of the time-stepping 
procedure to converge the unsteady solution at each time step. At each subiteration level, the 
discrete equations are approximately factored and solved using a modified two-pass matrix 
solver [Whi90], based on the Gauss-Seidel iteration procedure. 

In the present study, we have extended and applied TURBO to predict unsteady flows 
through single vibrating blade rows operating within cylindrical annular ducts. In particular, 
we have implemented a blade vibration capability into TURBO so that unsteady solutions 
can be determined over a domain that deforms with a vibratory blade rotation. In addition, 
we have incorporated far-field conditions into TURBO %o render computational inlet and 
exit boundaries transparent to outgoing disturbances and to allow incoming aerodynamic 
disturbances to be prescribed as approximate solutions to the governing equations. 

The goals of the present effort have been to extend, demonstrate and validate the TURBO 
analysis for blade flutter applications. In the present version of TURBO, the implicit, wave- 
split, finite- volume analysis, developed at MSU, is applied to predict the unsteady flow in the 
near field, and coupled at the computational inflow and outflow boundaries, to far-field eige- 
nanalyses for the unsteady perturbations of fully-developed, axisymmetric, mean flows. The 
resulting analysis is described in this report and applied to predict unsteady flows through 
a three-dimensional version of the 10th Standard Cascade [FV93] and the NASA Rotor 67 
fan. We have considered inviscid unsteady subsonic and transonic flows excited by prescribed 
blade vibrations and, for validation purposes, we have compared the TURBO results for 10th 
Standard Cascade with those based on the two-dimensional, potential-based linearization, 
LINFLO [Ver93], and those based on the three-dimensional, linearized Euler analysis, LIN- 
FLUX [MV97]. The predictions indicate that the current version of the 3D TURBO analysis 
can provide useful and accurate unsteady aerodynamic response information, provided that 
meshes of sufficient density and clustering are employed. 
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2. Unsteady Flow through a Vibrating Blade Row 


We consider the flow, at high Reynolds number (Re) and with negligible body forces, of 
a perfect gas with constant specific heats through a rotating and vibrating blade row that 
operates within a stationary annular duct (see Figure 1). The duct is of infinite axial extent 
and has hub and duct radii, r = r H and r = r D , respectively, and the blade row consists 
of N b blades which rotate about the duct axis at constant angular velocity fi = fie*. We 
assume that the flows far upstream and far downstream from the blade row are at most 
small perturbations of fully-developed, axisymmetric, steady background flows, and that in 
the absence of a vibratory motion, the blades are identical in shape, equally spaced around 
the rotor, and identical in orientation relative to the axisymmetric inlet flow. 

We will examine this unsteady flow in both stationary and rotating frames of reference, 
in terms of cylindrical (£, r, 9, t) and Cartesian (xi,X 2 ,X 3 , i) = (£, r cos 6, r sin 6, t) co- 
ordinates. Here £ and r measure distance along and radially outward from the duct axis, 
respectively, and 6 measures angular distance in the direction opposite to the direction of 
rotation, which is counterclockwise relative to an observer looking in the axial flow direction. 
When necessary, we will use the superscripts abs or rel to indicate that a physical quantity 
is measured relative to stationary or the rotating frame; e.g., 0 abs = 0 rel + fit. 

We intend to numerically resolve the unsteady flow, in terms of curvilinear spatial co- 
ordinates, on a computational grid that rotates with the blade row and deforms with the 
vibratory blade motion. The vector 72. (x, f) describes the displacement of a moving field 
(grid) point, x, relative to its reference or mean position, x, in the rotating frame. The 
displacement field, 72., is prescribed so that the solution domain deforms with the vibratory 
motions of the blades and is rigid far from the blade row. 

For aeroelastic and aeroacoustic applications, we are usually interested in a restricted class 
of unsteady flows; those in which the unsteady fluctuations can be regarded as disturbances 
to a background flow that is steady in a blade-fixed, rotating reference frame. Moreover, 
the steady background flows far upstream (say £ < £-) and far downstream (£ > £+) 

from the blade row can be assumed to consist of at most a small steady perturbation from a 
fully-developed, axisymmetric, steady flow. The time-dependent or unsteady fluctuations in 
these flows arise from temporally and circumferentially periodic unsteady excitations, i.e., 
prescribed vibratory blade motions and prescribed aerodynamic disturbances at inlet and 
exit that carry energy towards the blade row. 

For example, if the blades vibrate at reduced frequency, u, as seen by an observer in the 
rotating frame, and at constant interblade phase angle, cr , we can write 

72. Bn (f, 9 + 2t m/JV B , £, t ) = T n Re{R B (f , 6 , 0 exp[i(u;t + na)]}, x on B . (2.1) 

Here, 72. Bn is the displacement of a point on the nth moving blade surface from its mean 
position in the rotating frame; T n is a rotation matrix, which rotates vectors through n 
passages; n = 0,1,2,..., Nb — 1 is a blade index; Re{ } denotes the real part of { Rb 
is the complex amplitude of the reference (n = 0) blade displacement; and B refers to the 
mean position of the reference blade. The interblade phase angle, a, is determined by the 
nodal diameter pattern of the vibratory blade motion, i.e., a = 2 -kNq/Nb, where |iV D |, 
the number of nodal diameters, is the integer count of the number of times a disturbance 
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pattern repeats around the wheel. The sign of iV D is deteimined by the direction of rotation 
of the disturbance pattern. If the vibratory disturbance pattern moves in the direction of 
blade rotation, i.e., the negative ^-direction, then N-o > 0. We should note that for an 
excitation of the form (2.1) the unsteady flow will be periodic over N P blade passages, where 
N P = N b /\N d \, N d # 0. 

The unsteady disturbances in the far upstream and far downstream regions are, in part, 
prescribed as a fluid dynamic excitation and, in part, depend upon the interaction between 
the fluid and the blading. Typically, an unsteady aerodynamic excitation is represented by a 
linear combination of fundamental disturbances that are harmonic in time at relative reduced 
frequency w, and in the circumferential direction at angular wave number fh = N D + mN B , 
where m is an interger. The frequency of such an excitation in the stationary or absolute 
frame is = ui — mQ, where the term — roQ accounts for the Doppler shift. In the 
present study, we will restrict our consideration to unsteady flows driven by prescribed blade 
motions; therefore, all external aerodynamic excitations axe set equal to zero. 
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3. Fluid Dynamic Equations 


In the present discussion, all physical variables are dimensionless. Lengths are scaled 
with respect to the reference length L*\ time with respect to the ratio L* jV* where V* is the 
reference flow speed; velocity with respect to V * ; density with respect to a reference density 
p*\ stress, and therefore, pressure, with respect to p*{V*) 2 ] and specific internal energy with 
respect to (V*) 2 . The superscript * refers to a dimensional reference value of a flow variable. 
The scalings for the remaining variables can be determined from the equations given below, 
which have the same general forms as their dimensional counterparts. The reference length, 
L*, is typically taken to be the blade chord at the reference radial location the reference 
fluid density and flow speed, to be the temporally- and circumferentially-averaged inlet 
density and relative flow speed at r* = r*^, respectively. 


3.1 Governing Equations 

The field equations that govern the unsteady flow are determined from the conservation 
laws for mass, momentum and energy, the thermodynamic relations for a perfect gas, and 
the constitutive relations for a Newtonian fluid. After ensemble averaging these equations 
and applying an algebraic turbulence model, we arrive at the following form of the Reynolds- 
averaged Navier-Stokes equations 

?L + £ A+e,)=8, 

in which a summation over repeated indices is implied. 

The state, U, flux, tj and G„ j = 1,2,3, and source term, S, vectors in equation (3.1), 

are given by 


_ P 

pV XI 

pV* 2 

pVx 3 


F = 

r 3 


pV Xj 

pV Xi Vxj+PSij 

~pV x ,-V Xl + P5 2j 

pV_x 3 %+P$. 3i 
p(Et + p ipyv X j 


, Gj = n I2Ij . 

— ^■X 3 Xj 

. d* Qxj 


S — fl p( 2 Vi 3 ■+■ ^2^2) , 

-^(2 v; 2 - 

p£l{x 2 ^c 2 ^ 3 ^: 3 ) . 

where p,V, E T = E + V 2 /2 and P = (7 - l)~p{E T - V 2 /2) are the fluid density, relative 
velocity, relative specific total internal energy, and pressure, respectively, and 7 is the fluid 
specific heat ratio of the fluid. The components of the viscous stress tensor, II, and, assuming 
Fourier’s law for the conduction of heat, those of the heat flux vector, Q, are given by 


n XiX j — PefiP^ 


~/ x . dv x . dV k . 

_*L I ZL _ 9 —R 

x. dxi dxic 13 
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and 


(3.4) 


- ~ 11 UJL 

Q Xi = -KrfPr Re . 

Here f = 7 (Et - V 2 /2) is the fluid temperature, and Re and Pr are the Reynolds number 
and Prandtl number, respectively, of the flow. The inviscid flux, Fj, and source term vectors 
in equation (3.1) can also be written as explicit functions of the state variables; i.e., 


Fj(U) 


U »1 


0 

u j+l u 2 /u x + P 5 ij 


0 

Uj+iUs/Ui + PS* 

, S(U,x) = 

FPU\X 2 •+■ 2S1L^4 

Uj+iUi/Ui + P 5 3j 


tfthxs - 2QU Z 

. Uj^fa+pyu! J 


. tf(U 3 x 2 + V 4x3) . 


(3.5) 


where P = ( 7 - 1){U 5 - + (% + Ufll 2]. 

The matrix equation (3.1) describes the unsteady flow at a moving field point, x, as 
seen by an observer fixed in a reference frame that rotates at constant velocity, £2. If 
we set £2 = 0, we recover an equation that describes the flow in a stationary frame of 
reference. Equation (3.1) describes the behavior of the ensemble- or Reynolds-averaged 
values of the time-dependent flow variables. The effects of random turbulent fluctuations 
have been accommodated by using the effective viscosity, = p + t, and the effective 
thermal conductivity, R^ — p + (Pr/Pr T )l, where Pr T is the turbulent Prandtl number, 
in the definitions of the viscous stress tensor and the heat flux vector, respectively. The 
molecular viscosity, p, and thermal conductivity, R = p, are related to the temperature 
using Sutherland’s Law, the eddy viscosity, t, is determined using the Baldwin Lomax [BL78] 
algebraic turbulence model, and we set Pr = 0.72 and Ptt = 0.9. 


Transformation to Curvilinear Coordinates 


It is convenient to solve the foregoing field equations in terms of body-fitted curvilinear 
spatial coordinates (ai, 0:2, ct$) and the time r = t, where the positive directions of ai, ci!2 
and £*3 coordinate curves generally point in the streamwise, the spanwise (hub-to-tip) and 
the pitchwise directions. After introducing the transformation (x, t) -4 [a(x, t), r] into 
equation (3.1), we find that 


db_ 

dr 


+ 


a 


d 

datj 


(Fj + Gj) 



(3.6) 


where 

U = J" 1 U , 



' dotj - 9 q , \ 

^ u+ ^ F ‘j 



§ = J" 1 S 


(3.7) 


and J is the Jacobian of the transformation (x, t) -4 (o,r). The Cartesian components of 
fl and Q are determined from the relations 


H ijij 


fi eS Re 1 


da k 8V Xt | da k 8V X] 
dxj doi k dxi da k 


-2 


da w dY Xk 

dx k da m 



(3.8) 
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and 


(3.9) 


_ ,da dT 

Q Xi = -KesPr Re • 

Boundary Conditions 

For turbomachinery applications the field equations (3.1) or (3.6) must be supplemented 
by boundary conditions at the blade surfaces and at the duct walls, periodicity conditions 
over Np blade passages; e.g., V(r, 9+2ttNp/Nb, 0 = T ^ p 'V(r,6, £), and far-field conditions 
at the computational inflow (£ = £-) and outflow (£ = £+) boundaries. Since transient 
unsteady aerodynamic behavior is usually not of interest, a precise knowledge of the initial 
state of the fluid is not required. No-slip conditions, i.e., 

V abs = f2 x r + 7l Bn for and V abs = 0 for r = r H ,r D , (3-10) 

where V abs = V rel 4- ft x r, apply at the blade moving surfaces, B n , and at the stationary 
duct walls, respectively. In addition, either the heat flux Q • or the temperature T 
must be prescribed at such surfaces. Temporally- and circumferentially-averaged values 
of the total temperature, the total pressure and the flow angle are specified as functions 
of radius at the computational inflow boundary, i.e., at £ = £_, and the temporally- and 
circumferentially-averaged pressure is specified at the outflow boundary, £ — £-h consistent 
with radial equilibrium. In general, the unsteady fluctuations at inlet and exit that carry 
energy towards the blade row must also be specified; those that carry energy away from the 
blade row must be determined as part of the nonlinear unsteady solution. 


3.2 High Reynolds Number Approximations 


Thin-Layer Equations 


For most flows of practical interest, the Reynolds number (Re) is sufficiently high so 
that the viscous effects are concentrated within thin layers that lie along the blade surfaces 
and the duct walls (boundary layers), and extend downstream from the blade trailing edges 
(wakes). Such flows can be described by approximate field equations, known as the thin- 
layer, Reynolds-averaged, Navier-Stokes equations, leading to a substantial reduction in the 
computational resources needed to determine viscous unsteady solutions. The thin-layer 
equations are derived from (3.6) by assuming that streamwise gradients of the viscous flux 

terms, i.e., dG\/da i, are small, and hence, can be neglected. In addition, in the £*2 and 03 
directions, normal second derivatives of the velocity components and the temperature are 
retained, but mixed second derivatives are regarded as negligible. 

The field equations resulting from the foregoing approximations have the form 


au 

dr 


la 


+ ^ + da. 


a -tl a £ tl 
0 g 2 +7 t-g 3 = s , 


dO!3 


(3.11) 


-TL z TL 

where the column vectors G 2 and G 3 are the thin-layer approximations to the viscous 

flux vectors G2 and G3, respectively. The boundary conditions to be used in conjunction 
with equation (3.11) are the same as those discussed above for the full viscous equations. 
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Inviscid Flow (Re — > oo) 


The field equations that govern the fluid motion in the inviscid limit (Re — > oo), i.e., the 


Euler equations, 


at 

dr 



(3.12) 


are obtained from equation (3.6) by setting 6, = 0, j = 1,2,3 or from equation (3.11) by 

setting G = 0, j = 2, 3. In principle, the inviscid field equations must be supplemented 
by jump conditions that apply at vortex-sheet wakes, W m , and at shocks, Shm- However, 
the usual practice is to solve the inviscid field equations over the entire fluid domain, thereby 
capturing dis continuous wake and shock phenomena. The inviscid flow is then determined 
as a solution of the Euler equations subject to flow tangency conditions, i.e., 


(V abs — fi x r — • n = 0 for x 6 B n and • n = 0 for r = th , tb (3.13) 

at the blade surfaces and the duct walls, respectively. The periodic and far-field conditions, 
used in the inviscid approximation, are the same as those indicated previously for Navier- 
Stokes simulations. 


3.3 Solution Strategy 

We require numerical solutions to the foregoing nonlinear unsteady boundary- value prob- 
lems over Np blade passages, where Np — Nb/\Nd\ if |A r c| ^ 0 and Np = 1 if No = 0, to 
predict the unsteady aerodynamic responses of a blade row to harmonic and circumferentially 
periodic, unsteady excitations. In the present study, we will seek such solutions for inviscid 
unsteady flows by matching a wave-split, finite-volume analysis for the unsteady flow in the 
near field, i.e., in the region £_ < f < f+, to approximate solutions for the unsteady pertur- 
bations of fully-developed, axisymmetric, mean flows in the regions far upstream (£ < £_) 
and far downstream (f > £+) of the blade row. Thus, we will solve the nonlinear unsteady 
equation (3.12) in the near field, a linearized form of this equation in the far-field, and match 
the near- and far-field solutions at the computational inflow and outflow boundaries. 

The displacement field 71 is assumed to vary harmonically with time, i.e., 7t(x, t) = 
Re{R(x) exp(iu;t)}. The complex-amplitude of this field, R(x), must be prescribed over 
the entire solution domain. In the present study, R is defined so that the solution domain 
deforms with the blade motion (i.e., R = Rb„ f° r 5c € B n )> slides along the hub and duct 
walls (Rn = 0 for f = r H , td), and remains rigid far from the blade row (R = 0 for £ > £*)■ 
In addition, R(x) is prescribed along one blade-to-blade periodic boundary, such that it is 
continuous at the blade leading and trailing edges and decays exponentially away from the 
blade row. At the other periodic boundary, R is set so as to satisfy the periodicity condition 

R(f, 6 + 2i:Np/Nb, 0 = T*„R(f , 0, 0 . (3.14) 

In the near field, R(x) is first determined along the hub and duct walls as solutions of 
Laplace’s equation, — 0, in two dimensions. It is then determined in the interior of 
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the computational domain as a solution of Laplace’s equation in three dimensions, subject 
to Dirichlet boundary conditions given above. For the unsteady excitations being considered 
herein, it is sufficient to solve for the foregoing linear boundary value problems over a single 
extended blade-passage region, since R(x) can be specified in the remaining passages using 
the phase-lagged periodicity condition 

R(f, e + 27m/N B , 0 = T n R{r, 0, f)exp {ina) . (3.15) 

Also, note, that for unsteady flows in which no blades vibrations occur, one would simply 
set R = 0. 

The near-field, finite-volume analysis, which, at present, is performed in the stationary 
frame in ter ms of absolute flow variables, is described in §4 of this report. The fax-field 
eigenanalyses, which are performed in the rotating frame in terms of relative flow variables 
are described in §5. These have been coupled and implemented into the TURBO code, 
which is demonstrated via the numerical results for inviscid unsteady flows, presented in 
§6. Although our numerical results pertain only to inviscid flows, we have included viscous 
equations in this and the following sections to provide a framework for future work. 
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4. Near-Field, Finite- Volume Analysis 


A flux split, finite-volume analysis for nonlinear, inviscid, unsteady flows has been de- 
veloped [Jan89, JW89, JW90, JHW92], and implemented into the turbomachinery unsteady 
flow code, TURBO. This analysis was later extended for the prediction of viscous flows 
[CW93, CCA94]. TURBO is an implicit, multi-block, cell-centered, finite- volume code, that 
can be used to predict three-dimensional, nonlinear, inviscid and viscous, steady and un- 
steady flows through and around blade rows. The fluid dynamic equations are solved in 
a stationary reference frame over a solution domain that rotates with the blade row and 
deforms with the vibratory blade motions. A brief description of the TURBO analysis is 
given below. Additional information can be found in the references cited above. 

The computational mesh used in TURBO is a sheared H-mesh. This structured mesh 
defines a curvilinear coordinate system, in which the coordinate curves lie along the bound- 
aries of the physical domain, such that there is a one-to-one correspondence between the 
points, x, in the physical domain and the points, a, in the computational domain. A time- 
dependent coordinate transformation, (x,t) — » (a,r), where x = x + 7Z.(x, t), from the 
rotating physical domain, in which the grid deforms with the blade motion, to a computa- 
tional domain, in which the grid is stationary, uniform, and orthogonal, is applied to simplify 
the implementation of numerical differencing and flow boundary conditions. The on, a? and 
0:3 computational coordinates, or the /, J, K computational mesh indices, refer to the axial, 
spanwise and pitchwise directions, respectively. Cell faces are surfaces of constant computa- 
tional coordinate, so that each cell is bounded by the six surfaces: c*i = I — 1/2 and I + 1/2, 
and c*2 = J — 1/2 and J + 1/2, and 0:3 = K — 1/2 and K 4- 1/2. 

4.1 Finite- Volume Equations 

For a finite- volume discretization of the governing field equations, the time-dependent 
geometrical properties of the mesh cells in physical space are required. These include the 
cell volume, ~d = J -1 , the volume swept out per unit time by the constant a y face as the cell 
interface moves, •dj = J~ l daj/dt, and the area of the constant ctj cell face projected in the 
Xk direction, Aj k = J~ l dotj / dx k . These geometric properties of a cell axe determined from 
the instantaneous locations of the cell vertices in physical space. 

The finite-volume spatial discretization [Jan89, CW93] of equation (3.5), expressed in 
the stationary frame (f2 = 0), can be written as 

d\J/dr = —6jFj - 6jGj = -R (4.1) 

where U = tfU, F, = — djU -(- Aj k F k and Gj = Aj k Gj Here, U represents an average of 
the physical state vector over a cell volume; Fj is the inviscid or convective flux and G,- is 
the viscous or diffusive flux, across a constant aj cell face; and R is the residual. The flux 
vectors F and G depend on the physical state varibles and the cell surface properties, and 

the residual R is a nonlinear function of the physical s tate vector, U. The operator Sj in 
equation (4.1) denotes the difference in the Oj-direction across adjacent cell interfaces, and 
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the repeated j index implies summation over all computational coordinate directions, so that 
the terms SjFj and 8j Gj are the net inviscid and viscous fluxes, respectively, through a cell. 

The time-derivative in (4.1) is approximated using a second-order, implicit, three-level, 
backward, difference approximation. After applying this scheme and separating the time 
dependence of the state vector and the cell volume, we find that 

i?AU n + r" + 1 = S n+1 = $ n-1 AU n_1 /2AT - U n (3$ n+1 - 4$ n + $ n_1 )/2A r , (4.2) 


where d = 3$ n+1 /2A r, AU n = U n+1 - U n , the vector Q n+1 is determined by the terms on 
the right-hand side of (4.2), and the superscript n refers to the nth time level. The nonlinear 
equation (4.2) is solved at each time step using a Newton iteration procedure in which the 
viscous flux terms are treated explicitly, i.e., we set 


d AU* + S s (dFj/dV i A U 1 ) = - 0") -R + S“ +i 


(4.3) 


where p = 1,2,... , is the Newton iteration index, U° = U n , AU P = U p — U p-1 _, and U p is 
the Newton update to the state vector. Once the Newton iteration converges AU P = 0, and 
U p = U n+1 . 


4.2 Evaluation of Flux Terms 


To simplify the description of the spatial discretizations that are used to evaluate the 
flux terms that appear in (4.3), we consider a “one-dimensional inviscid flow” in which, for 

example, F j = F is the inviscid flux vector in the aj = a computational coordinate direction. 
The subscript J will refer to the cell volume bounded by the cell surfaces at a = J + 1/2 
and a = J - 1/2. Extensions of the equations that follow to multi-dimensional flows are 
straightforward conceptually, but involve the use of tedious additional nomenclature. 

A cell-centered finite-volume discretization requires that flux information at say the J+l 
cell interface be computed in terms of the values of the state variables, Uj and Uj +1 , in the 
neighboring cell volumes and the geometric properties of the grid, i.e., $ 7 + 1/2 and A J+ 1 / 2 , at 
the cell interface. In the TURBO analysis, a flux splitting technique is applied to evaluate 
the interfacial inviscid fluxes. It is based on a similarity transformation and an eigenvalue 
decomposition, of the flux Jacobian matrix, dF/dU, into matrices that account for right (+) 
and left (— ) traveling disturbances. 

Thus, the matrix 3F/3 U|u j Gj+i/2 , where the subscripts indicate that the flux Jacobian 
matrix is evaluated in terms of the state vector in Jth cell and the surface metrics, indicated 
by Gj+ 1 / 2 , at the (J + l/2)th cell interface, is split according to 


dF 

_ dF 

+ ^ 

dF 

+ “77 

= f(A + + A )T 1 

au 

~ au 

Uj,C?j+i/2 

d\J 

Uj,Gj + i/2 

Uj»Gj + 1 /2 


Uj,Gj + i/2 


(4.4) 


where the matrices T and T contain the right and left eigenvectors, respectively, of 9F/5U, 
and A + and A are diagonal matrices containing the positive (+) and negative (-) eigenval- 
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ues. The eigenvalues of the flux Jacobian matrix are used to determine which characteristic 
modes are taken into account, thus controlling the direction of spatial differencing. 

Two different approximations are applied to evaluate interfacial, inviscid, flux informa- 
tion. One is based on flux vector splitting [SW81], and is used to evaluate the left-hand 
side flux ter ms in (4.3). This approximation is only first-order accurate, but allows for a 
convenient approximate factorization of equation (4.3), which facilitates the time-marching 
solution. The other is based on flux difference splitting [Roe81], and is used to evaluate the 
right-hand side inviscid flux terms. In TURBO, Roe’s first-order accurate, flux-difference 
splitting approximation is modified, by adding corrective fluxes, to achieve higher order 
spatial accuracy. 

Left- Hand-Side Flux Terms 


The flux- vector splitting approximation to the Newton update to the inviscid flux vector 
at the J + 1/2 cell interface is 





AU P 

up-> / 


J+l/2 


dF 

au 


AU$ + 


*r l .o 5 +i/» 


dF 

dV 


AU 5 +1 , 


(4.5) 


TT P “ 1 

u J+l’ L, J+l/2 


where the subscripts on the right-hand-side, flux Jacobian matrices indicate that these ma- 
trices are evaluated in terms of the state vector, U p— 1 , in the indicated cell volume, J or 
J + 1, and the swept volume and surface area at the n + 1 time level and the J+l/2 cell 
interface. The approximation (4.5) results in first order spatial accuracy, but it is only used 
to construct an approximate factorization of the Newton iteration equation (4.3). Therefore, 
any errors that are introduced, do not appear in the converged final solution. 

The terms in (4.5) are spatially differenced in the a-direction to determine the Newton 
update to the net inviscid flux through the Jth cell volume; i.e, 



dF 

au 


VVp-1 grn + 1 
U J ’ U J + 1 


3F 

AU P + ~ 

J au 


J+l/2 


au; +1 


TJ P ~ 1 <7 n+1 
U J+l’ U J+l/2 


aF 

dU 


1 + 


tjp- i G 71 * 1 

u J-l > U J-1 


dF 

AU^ - -- 

J 1 au 


1/2 


au; , 


ttP ” 1 G n+1 
U J >Ur J- 1/2 


(4.6) 


to determine the net flux through the Jth cell volume. 


Right-Hand Side Inviscid Flux Terms 

The inviscid flux vectors that appear in the residual on the right-hand side of (4.3) are 

evaluated using flux difference splitting [Roe81]. In TURBO, the flux, Fj+i/ 2 ; at the J + 1/2 
cell interface is evaluated in terms of the flux in the cell to the left ( J) of the interface and 
the flux due to waves approaching the interface from the right. Thus, we set 


F J+l/2 = F(Vj,Gj +1/2 ) + 


aF 

au 


(Uj+i - Uj) , 




(4.7) 
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where F(Uj,Gj + i/ 2 ) is a flux based on the state vector in the Jth cell and the cell metrics, 
it and A, at the J + 1/2 cell interface, and the flux Jacobian matrix 9F/9 U|ur« /2 Gj+1/2 

is evaluated in terms of the intermediate state vector, U , and the cell metrics at the 
J + 1/2 interface. The intermediate state vector, U^. e 1/2 , is defined using the relations: 


and 


fiZz, = \!hh « . v «?', /2 = 


\/a/Vj + y/p j + t V j+l 


17 Roe 
^TJ+l/Z 


Pj + \Jpj + i 

yffhEr,J + 


(4.8) 


The discrete approximation (4.7) is first-order accurate, since the interfacial fluxes are 
based only upon information from adjacent cells. Higher order spatial accuracy can be 
achieved by adding corrective fluxes to the right-hand side of (4.7), which bring in information 
from additional neighboring cells. The corrective perturbation flux at the J + 1/2 interface 
is comprised of right traveling waves at the upstream interface, J 1/2, of the </th cell and 
left traveling waves at the downstream interface, J + 3/2, of the ( J + l)th cell. These waves 
are approximated using the Roe-averaged, flux Jacobian matrix at the J + 1/2 interface. 
Thus, the enhanced approximation to the perturbation flux is obtained by adding terms of 
the form 


1 8F_ 

2 d\5 


+ 


U ^* a / 2-^+ 1 / 2 


(Uj - Uj-i) 


and 


1 5F 

2 au 


(U J+2 - Uj+i) 

<% /2 , *7+1/2 


to the right-hand side of (4.7), which should result in second order spatial accuracy. Flux 
limiters [VL74] are used in conjunction with the corrective fluxes to control dispersive errors, 
such as those that occur at shocks and at stagnation points. The limiters are activated by 
changes in sign in the jumps in the characteristic variables at adjacent interfaces. 

Once the interfacial fluxes have been computed, they are spatially differenced, i.e., 


8 F 


= F 


.7+1/2 


-F, 


J- 1/2 


(4.9) 


to compute the net inviscid flux through the Jth control volume. A second-order discrete 
approximation is used to evaluate the interfacial fluxes in (4.9). Note that in computing the 
residual at the p - 1 Newton iteration level, the flux vectors in (4.9) are evaluated in terms 
of the state variables at the p — 1 iteration level and the grid properties at the n + 1 time 
level. 


Right-Hand Side Viscous Flux Terms 

At each step of the Newton iteration procedure, the viscous flux vector, G, is evaluated 
at the cell interfaces in terms of the values of the flow variables, at the p - 1 iteration level, in 
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the cell volumes to the left and right of the interface and the area of the cell interface at the 
n + 1 time level. The individual terms that make up this vector are considered separately. 
Derivatives of the fluid properties at cell interfaces are evaluated in terms of the property 
values in the cell volumes to the left and right of the interface, using central difference 
approximations; the velocities at cell interfaces, by averaging their values in the adjacent 
cell volumes. Once the viscous fluxes at the cell interfaces have been computed, they are 
spatially differenced according to 


6G 


j 


— Gj+i/2 — Gj_i/2 


(4.10) 


to compute the net viscous flux through the cell volume. 


4.3 Solution Procedure 


The spatial difference approximation (4.6) leads to an approximate factorization of the 
Newton iteration equation (4.3) of the form 

fjjAu; - mLau;., + jviJ +i au5 +1 = -h(V7' - uj) - r 7‘ + 0" +1 , (4-n) 


where the index J refers to the «7th computational cel' and the D and M matrices are 
evaluated based on the state vector U p_1 . The D matrix contains the diagonal elements of 

z 4 - z — 

the iteration matrix, and the M and M matrices contain the off-diagonal elements in the 
negative and positive computational coordinate directions, respectively, i.e., 


Dj 


2 + 


Mj-i 


1 + 


dF 

dV 


+ 


fjp-i /7 n+1 
u J ’ u jr+ 1/2 


dF 

dV 




1/2 


dF 

9U 




dF 

and M , , , = — ^ 
J+1 dU 


fjP-1 G n+1 
u ;+ p u ;+ i /2 


(4.12) 


To reduce the errors introduced by the approximate factorization, equation (4.11) is 
solved for AU P using a symmetric Gauss-Seidel subiteration procedure. The first subit- 
eration is over positive grid indices; the second, over negative grid indices. The subitera- 
tion procedure is thus an LU decomposition of the Newton iteration matrix, with forward 
and backward substitution. Once the Gauss-Seidel subiteration procedure converges, equa- 
tion (4.3) is satisfied, and the calculation proceeds to the next Newton-iteration level. As the 
solution at time r = r n+1 converges, any errors introduced by the Newton iteration or the 
approximate factorization vanish. Only the errors in the calculation of the residual of equa- 
tion (4.1) remain. The terms that make up this residua are calculated using second-order 
accurate difference approximations. 
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Boundary Conditions 


The field equation (4.1) must be solved subject to conditions at the boundaries of the 
near-field computational domain. The flow tangency conditions used in the inviscid version 
of TURBO, cf. (3.13), are implemented by using phantom cells inside a solid surface. The 
density and pressure in a phantom cell are defined by a first-order accurate reflection con- 
dition, and the phantom cell velocity is defined such that the velocity at a solid surface, 
which is the average of the velocities in the phantom and the interior cells, satisfies the flow 
tangency condition, in a manner consistent with the finite volume discretization. Periodicity 
conditions; e.g., Vy = T^ p "V" l , where the subscripts U and L refer to the upper and lower 
periodic boundaries of the computational domain, are imposed at the pitchwise boundaries. 
Finally, as discussed in the next section, analytic/numeric far-field solutions, based on re- 
duced forms of the governing equations, are matched to the numerical near-field solution at 
the computational inflow and outflow boundaries (£ = £t) - 

The current TURBO implementation uses explicit boundary conditions, which are incor- 
porated into the SGS iteration procedure, so that the boundary conditions are imposed in a 
semi-implicit manner. This treatment has been found to yield better convergence properties 
than a purely explicit implementation. 
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5. Far-Field Eigenanalyses 


Fax-field solutions, based on reduced sets of governing equations, can be applied to restrict 
axial extent of the near-field computational domain. To develop such solutions, we require 
an inviscid form of the field equation (3.1), that applies at fixed locations (x = x) in the 
rotating frame. Expressed in terms of fluid dynamic variables, V and E T , measured relative 
to a reference frame that rotates with the blade row and in terms of rotating cylindrical 
coordinates, this equation has the form 


0U .drtr x dF 6 0F € _5 

In +T ~aT + r ar + IT 


where the state and source-term vectors in equation (5.1) are given by 


f ~p 1 


0 


0 

PVr 


{Uzf/U.+P 


2 U 3 + UiQr 

~pV e 

(/)l 

It 

1 

i-» 

-UlUz/ih 

> + Q. < 

-2 U 2 



0 


0 

. P-E’r > 


0 


k U&r , 


(5.1) 


(5.2) 


The flux vectors F r (U), F<?(U) and F f (U) and the pressure P(U)_have functional forms 
similar to those indicated previously for the Fj-(U), j = 1,2,3, and P(U) in §3.1. 


5.1 Unsteady Perturbations in the Far Field 

To determine approximate solutions to equation (5.1), that describe the flows far up- 
stream (£ < £_) and far downstream (£ > £+) of a blade row, we first expand the unsteady 
state vector, U, into an asymptotic series of the form 

U[x, t] = U(x) + u(x, <) + ...= U(x) + Pe{u(x) exp(iurf)} + . . . , (5.3) 

where the column vectors U(x) and u(x, t) contain the conservation variables for the zeroth- 
order background flow, which is steady in the rotating frame, and the first-order unsteady 
perturbation, respectively, and the dots refer to higher ol der terms. The components of the 
vector u are the complex amplitudes of the first-order unsteady conservation variables, i.e., 
u T = [p, pv r + pV T , pve+pVg , pv^ + pV^, pe T +pE T ] where p, V and E T and p, v, and e T are 
the steady and the complex amplitudes of the first-order unsteady, primitive, flow variables, 
respectively. The unsteady flux, say F r , and source term, S, vectors are approximated using 
Taylor series expansions about the mean flow state, U, i.e., 

F,(U) = F r (U) + §u + . . . and S(U) = S(U) + + . . . . (5.4) 

Field equations, that describe the steady and the first-order unsteady flows in the far 
upstream and far downstream regions, are determined by substituting the foregoing series 
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expansions into the nonlinear, time-dependent equation (5.1), and equating terms of like 
order. The resulting equations for the zeroth- and first-order flows are 


_ : drF r , ^dFe , dF * 
hr — H 


dr 


dO di 


= S 


(5.5) 


and 


,d(rAu) _! 9Bu dCu 

!wu+r -V i+r V + ir ~ Du=0 


(5.6) 


respectively, where A = 3F r /dV, B = 3F S JdV and C = 8F ( /3U are flux Jacobian 
matrices and D = dS /dU is the source-term Jacobian. 

We assume that, far from the blade row, the mean or steady flow quantities are dependent 
only on radial position; i.e., p = p(r ) , P — P{r), etc., and that the radial component of the 
steady velocity is negligible; i.e., V = Vg (r )e$ + V^(r)e^. Under these conditions, the steady 
field equation (5.5) reduces to 


p- 1 ^- = r ~ ly e + mVe + ft 2 r = r" 1 ^ 8 ) 2 , (5.7) 

dr 

where y 0 abs = V$ + x r is the absolute circumferential velocity. We also assume that the 
kinematic and thermodynamic data needed, in conjunction with (5.7), to completely specify 
the steady background flow in the far field axe available. 

For the mean flow conditions just described, the linearized unsteady equation (5.6) re- 
duces to .. . V Q O 

iuu + r- 1 +r- 1 B 3+C 2 ^-Du = Q, (5.8) 

or ov a £ 

where the subscript 2 on the Jacobian matrices in (5.8) indicates that they are evaluated at 
U 2 = pV T = 0; e.g., A 2 = dF r /dV\u,= 0 . 


Uniform Mean Flow 

For the special case of a uniform mean flow, in the absolute frame, i.e., V£ hs = 0 and 
is a constant, an exact solution can be determined for the first-order unsteady perturbation. 
This solution indicates that an arbitrary unsteady disturbance can be represented as the 
sum of independent entropic, vortical and irrotational acoustic disturbances. A state vector, 
representing a linear combination of entropic and vortical disturbances is a solution 
of (5.8) that satisfies the convection equation Duc/Dt = ( iui -I- V • V)uc = 0. Thus, such 
disturbances are convected by the mean flow, and, for an unsteady flow occurring at temporal 
frequency u in the rotating frame, u c has a general solution of the form 

u c = J2 u™(r) exp[i(K C>m £ + m0)] . (5.9) 

m=-oo 

Here, the are arbitrary functions of radius, fh = N-^+mNs and = -(w-mfi x )V^ 1 = 

are t h e circumferential angular and axial linear wave numbers of the mth distur- 
bance, and u4 bs = u ; - mfi x is the temporal frequency of the mth disturbance as seen by an 
observer fixed in the absolute frame. 
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The acoustic disturbances are governed by a convected wave equation for the unsteady 
pressure, which can be solved analytically [TS62]. The resulting solution for the complex- 
amplitude of the unsteady pressure in a subsonic axial mean flow is 

00 oo 

P= Y e ™ e £ [*w eX P(*m^) + Pm* ex P(Xm#i£)] ( r ) • ( 5 - 10 ) 

m=— oo jx=0 

Here p* are the amplitudes of the downstream and upstream traveling pressure waves, and 
E m(t {r) = Jfhikmtxr) + Qm^Ymikm^r) are the “characteristic E-functions” of [TS62]. The 
E-functions are combinations of Bessel functions, of order fn, of the first and second kinds. 
The constants k mfi and Q mt i in (5.10) are determined by the duct-wall boundary conditions, 
and the index p = 0, 1,2, ... indicates the number of zero crossings or nodes in the pth 
radial mode. 

The axial exponential coefficients, Xmn — Pm* + are given by 

xl „ = (l - Mf)-' A- 1 T ^(1 - Ml)AV w - ] , (5.11) 

where = V^/A < 1 is the mean axial Mach number. If (w^ 3 ) 2 > (1 - Ml)A 2 k then 
the Xmp 3X6 purely imaginary, and the mp th pressure patterns propagate. If the Xmn 3X6 
complex, then one pattern attenuates, and the other grows exponentially, with axial distance. 
The perturbation state vector for acoustic disturbances is given by 

n A = Y e ifh0 Y [% exp + u^exp (xm^)] > ( 5 - 12 ) 

m=— oo fi=0 

where the modal state vectors, u^ M (r), are determined, in terms of the pressure, from the 
linearized unsteady field equations. Note that, in addition to different axial behaviors, the 
acoustic disturbances in (5.12) and the convected disturbances in (5.9) have different radial 
behaviors. The former occur in radial modes, the shapes of which are determined by the 
unsteady field equations, whereas the latter have arbitrary radial dependence. 

Nonuniform Mean Flow 

Guided by the exact solutions for uniform mean flows, approximate solutions to the 
linearized unsteady equation (5.8) can be constructed for nonuniform steady background 
flows [MV97]. For this purpose, we set u — + uw, where u c describes the convected 

disturbance field, and 

OO OO 

u w = Y, exp(im0) Y u£ n (r ) exp(* mn f) , (5.13) 

m=-oo n=0 

describes a series of modal type disturbances. The summation over n in (5.13) includes 
different possible types of modes, such as upstream or downstream traveling modes as well 
as modes with different numbers of radial zero crossings. 

The convected disturbance is a solution to equation (5.8) of the form (5.9), but, for 
nonuniform mean flows, the axial wave numbers, ^ >m (r) = — [oj + fhr~ l Ve(r)]/V^(r) , are 
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functions of radius. The convected field may contain entropic and/or vortical disturbances, 
depending on the properties of the mean flow. However, for a general nonisentropic, rota- 
tional mean flow, no convected field will exist. 

The modal disturbances are determined by substituting the assumed form of the solution, 

(5.13), into the field equation (5.8), yielding the system 

iul u£ n + r" 1 ^ (r A 2 u£ n ) + iffir' 1 B 2 < n + XmnC 2 u£„ - D 2 u£ n = 0 , (5.14) 

which must be solved numerically. After discretizing (5.14), and applying the hub- and 
duct- wall boundary conditions, v T = 0 for r = rn,rD, we obtain the matrix equation 

(P - XmnC 2 ) ll£ n = 0 (5-15) 

where P = -iwI-L(r, A 2 ) -zmr' 1 B 2 +D 2 and L(r, A 2 ) is a finite difference approximation 
to r -1 c)(rA 2 uD/9r. The column vector ul in equation (5.15) contains an entry for each 
of the five conservation variables at each radial discretization point. 

Equation (5.15) can be solved [MV97] using a standard linear algebra routine, to de- 
termine the axial eigenvalues, Xmn, and the right eigenvectors, u£ n (r), of the modal far- 
field unsteady perturbations. The left eigenvectors are determined by solving the equation 
/p _ v C) H uL = 0, where the superscript H denotes the conjugate transpose. An or- 

thonormal set of left eigenvectors is obtained by setting (v mn ) = (Umn) Wll^nn/ '-' u mnJ- 

By invoking the orthogonality of left and right eigenvectors, the complex amplitudes of the 
modal disturbances, o mn , are determined by taking inner products involving v mn and Uw, 

Omn = <vl, ^ J* +2 * /NB Uw exp [ (Xmn£ + ifhff))dff ) . (5.16) 

5.2 Classification of Unsteady Disturbances 

Unsteady perturbations of uniform mean flows can be represented as superpositions of 
convected entropic and vortical disturbances and upstream- and downstream traveling irro- 
tational pressure disturbances. For nonuniform mean flows, the situation is more compli- 
cated [Kou95] . In particular, for rotational, but isentropic, mean flows, the unsteady entropy 
is an independent convected disturbance. However, because of the coupling between vortical 
and acoustic disturbances, due to mean-flow vorticity, neither convected vortical nor irrota- 
tional acoustic disturbances exist. Instead, nearly-convected or vorticity-dominated modal 
disturbances, that contain pressure, and acoustic or pressure-dominated disturbances, that 
contain vorticity, occur [GA96]. These types of disturbances emerge as solutions of the 
eigenvalue problem (5.14). 

The group velocity, 

v 9 , mn = du/dXmn = {yin, Cu£ n )/(v£ n , (dP/du) u£ n > , (5.17) 

i.e., the axial velocity at which an mnth modal disturbance carries energy, is used to classify 
modal disturbances. Nearly-convected disturbances travel downstream, without attenuation, 
at axial speeds slightly less than and slightly greater than the mean flow speed. If the mean 
axial Mach number is subsonic, acoustic disturbances travel both upstream and downstream. 
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For nonuniform mean flows, we can further decompose the unsteady state vector to 
account for the two types of modal disturbances. Thus, we set uiv = u^ + uat, where u a and 
Uat are the complex amplitudes the acoustic and the nearly-convected unsteady disturbances, 
respectively. The state vector for the acoustic disturbances has the form 

00 OO Pi 

iu(r, 0, 0 = X exp(im0) X «*p(x^aO + a m^,A^,A expfoj^*)] > 

m=— oo 

(5.18) 

where n indicates the number of radial nodes, and the — and 4- superscripts refer to down- 
stream and upstream traveling disturbances. The nearly-convected disturbances, i.e., 

oc °o n 1 

ujv(r, 0, 0 = X exp(im0) X [ ( W,w u m^( r ) «P(**£.nmJvO + exp(i/cJ m#t>Ar Oj • 

m=-oo ji=l 

(5.19) 

are also ordered by the number of radial nodes, but in this case starting with fj, = 1, and 
the — and + superscripts in (5.19) refer to disturbances that travel downstream at speeds 
slightly slower and slightly faster than the convection speed. 

In n um erical calculations, the series in equations, (5.18) and (5.19) must be truncated, 
since only finite number of circumferential and radial modes can be accommodated. Also, 
the numerical solutions to (5.15) will yield spurious modes; i.e., modes that satisfy the 
difference equation (5.15) but not the differential equation (5.14). The spurious modes must 
be elimina ted or filtered out, to yield a valid solution set. The filtering is based on the 
number of radial zero crossings or nodes and the point-to-point oscillations of a computed 
radial mode. Such criteria have been usually found to yield only genuine modes, but the 
filtering algorithm is still under development. As another caveat, since only a finite number 
of modes are retained after the truncation and filtering processes, the numerical far-field 
modal description may be incomplete. Based on previous work [MV97], the inclusion of 
spurious modes, or the exclusion of genuine modes, can be detrimental to both the accuracy 
and convergence of numerical solutions. 

5.3 Near-Field/Far-Field Matching Procedure 

The far field solutions must be applied in conjunction with a numerical near-field so- 
lution to determine the unsteady flow. Incoming disturbances (excitations) are prescribed, 
and outgoing disturbances are determined by matching the near- and fax-field solutions. 
The amplitudes of the outgoing modal disturbances are determined by taking inner prod- 
ucts, cf. (5.16), using the near-field state vector, u, in lieu of u w . This requires invoking 
the assumption that (v£ n ,u) « (v4 )n ,uw), i.e., that the left eigenmodes of the modal 
disturbances are orthogonal to the convected disturbances. Once the amplitudes of the out- 
going modes are determined, the wave- type modes axe superposed to provide a solution for 
uw = u a + u ff for a finite number of modes. 

At the upstream far-field boundary, the convected disturbance is set to describe any 
incident convected gust. At the downstream boundary, the wave-type modes are subtracted 
from the total unsteady disturbance and the remainder, uc = u — Uw, is treated as a 
convected disturbance. The convected disturbance in the region f > £+ is computed, by the 
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method of characteristics, as a solution of Du c /Dt = 0. Since the mean radial velocity has 
been assumed to be negligible, u c (r, 9, £) = u c {r, 9, f+) exp[-iu>(f - C+)/Vj] a 1011 ? constant- 
radius characteristics. 

In the near-field, the nonlinear unsteady equations are solved using the time-marching 
technique described in §4. After a pre-determined number of time steps, say N T , of the 
near-field solution, the amplitudes of the wave- type modes, i.e., a^ A and a^ jjV , and the 
complex amplitude of the far-downstream convected disturbance, uc(r,9,£+), are updated. 
The far-field solutions, which are the sums of wave-type and convected disturbances, are 
then updated, and used to supply the far-field boundary information needed for the next set 
of Nt near-field time-steps. 


23 



6. Numerical Results 


Unsteady aerodynamic response predictions will be presented to demonstrate the current 
capabilities of the TURBO code. First, we will consider subsonic unsteady flows through 
a rotor, based on the Tenth Standard Cascade Configuration [FV93], which is referred to 
in [MV97] as the 3D 10th Standard Cascade. Second, we will analyze the NASA Rotor 
67 [SWHS89], which is a research transonic fan consisting of 22 blades. 

We consider unsteady flows that are excited by prescribed single-degree-of-freedom, har- 
monic, blade motions (e.g., see Figure 2). The motions to be considered are pure translations 
normal to the sectional blade chords (R B = hn), and pure rotations about axes at the blade 
midchords (R B (x B ) = a x (x B - x P )). The complex amplitudes, h and a, of the bending 
and torsional vibrations are assumed to be constant along the span; n(r) = n g eg + is 
the unit normal to the blade chord at radius r, which is tangent to the cylinder r = constant; 
and x B -x P is the distance, at constant radius, to the point, x B (r), on the mean or reference 
blade surface from the point, x P (r), at the mean position of the torsional axis. 

The blade motions are termed subresonant if all fundamental acoustic response distur- 
bances attenuate with increasing axial distance from the blade row; superresonant (m, /i) 
if m and ji such disturbances persist in the far upstream and far downstream flow regions, 
respectively, and carry energy away from the blade row; and resonant if at least one acoustic 
response disturbance persists in either the far upstream or far downstream regions of the 
flow and carries energy along the blade row [Ver89b]. 

The TURBO analysis has been applied to predict unsteady surface pressure and the local 
(wc) and global (Wc) work per cycle responses to the prescribed blade vibrations. The local 
and global works per cycle are determined from the relat ions 

wc(x-b) = -v~ l f* P B ^^- • n B d(ut) and W c = (j^wc{x-B)dA B . (6.1) 

In equation (6.1), P B is the pressure acting at the point x B on the moving reference blade 
surface B, 71b is the displacement of this point relative to its mean position in the rotating 
frame, n B is a unit vector normal to B and pointing into the fluid, and dA B is a differential 
element of surface area. 

In addition to the nonlinear TURBO results, for purposes of comparison, we will also 
present response predictions for the 3D 10th Standard cascade based on the two-dimensional 
linearized analysis, LINFLO [Ver93] and three-dimensional linearized analysis, LINFLUX 
[MV97]. In LINFLO, the unsteady flow is regarded as a Email perturbation of a nonuniform, 
potential, steady background flow. The full-potential analysis CASPOF [Cas83] has been 
used to provide the steady background flows for the LINF LO calculations. In LINFLUX, the 
unsteady flow is regarded as a small perturbation of a nonuniform, Euler, steady background 
flow. The TURBO analysis has been used to provide the steady background flows for the 
LINFLUX calculations. 

The TURBO nonlinear steady-state solutions are determined, over a single extended 
blade passage, on an H-type grid. Because of our assumptions regarding the flows far from the 
blade row, the axial extent of the mesh must be chosen such that the mean flows at inlet and 
exit are at most small perturbations from steady background flows that are fully-developed 
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and axisymmetric. We have retained first-order steady quantities in our far-field expansions 
of the meanflow quantities to accommodate small axial and circumferential variations in 
the steady flow. This allows some flexibility in restricting the axial extent of the near-field 
computational domain. For the numerical examples presented in this report, we have found 
that an extent of one axial chord both upstream and downstream of the blade row to be 
conservatively sufficient. 

Since the TURBO code is written in terms of absolute-frame variables, the steady-state 
solutions for both numerical examples are obtained by marching the calculations in a time- 
accurate manner. These steady-state solutions are then used as inputs for both the unsteady 
TURBO computations as well as the linearized LINFLUX computations. The unsteady 
TURBO solutions are computed over single or multiple blade passages, depending on the 
interblade phase angle. The H-grids used for the present TURBO and LINFLUX calculations 
have been generated using either the IGB [BH92] or the TIGER [SS91] grid generation 
packages. 

For the 3D 10th Standard Cascade, the steady background flow at inlet is axial and 
uniform relative to a space-fixed or inertial reference frame. Thus, the absolute inlet Mach 
number, is a constant. For the Rotor 67 fan, the relative inlet Mach number 

is supersonic near the tip and subsonic near the hub. Thus, the steady background flow at 
inlet is only approximately axial and uniform relative to a spaced-fixed reference frame. 


6.1 3D 10th Standard Configuration 


The 3D Tenth Standard Cascade consists of 24 blades, which are twisted to reduce the 
variation in mean incidence due to blade rotation. The blades rotate within a cylindri- 
cal annular duct of inner radius r H = 3.395 and outer radius r D = 4.244. At midspan 
(r = r mid ), the blades are staggered at ©(r mid ) = 45 deg with a circumferential spacing, 
G(r„ii d ) = 27rr rnid /A r B , of unity, and the midspan blade section is a NACA 5506 airfoil, al- 
tered slightly [Ver89a] to have wedge-shaped trailing edges. The blade mean chord lines are 
located at 

rd = £tan© + nG(r), 0 < £ < cos©, n = 0, . . . , AT B — 1 , (6.2) 


where 


tan©(r) _ r 
tan©(r m i d ) r m j d 


(6.3) 


The axial chord is constant, hence, the leading and trailing edge £ and 9 coordinates are 
constant along the entire span. The airfoil chord varies from 0.946 at the hub to 1.057 at 
the tip, because of twist, and the local thickness to chord ratio varies to maintain constant 
thickness. The cascade operates in a uniform axial inlet flow, which occurs at M** = 0.4015, 
and rotates at an angular speed of |D| = 0.2145. This 3D configuration was chosen to match 
the subsonic Tenth Standard Configuration [FV93] at midspan, where the relative inlet Mach 
number, M_ is 0.7 and the relative inlet flow angle, Q-oo, is 55 deg. 

The H-grid for the 3D Tenth Standard Cascade consists of 141 axial, 41 tangential and 
11 radial surfaces (56,000 cells), and extends one axial chord upstream and downstream from 
the blade row. This is identical to the grid used in the 3D LINFLUX studies of [MV97]. 
Axial grid points are clustered near blade leading and trailing edges; circumferential grid 
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points, near blade surfaces; and radial points are distributed uniformly. In particular, the 
normal and chordwise grid spacings at a blade leading edge are 0.02% and 0.10% of chord, 
respectively, see Figure 3. 

The axial extent of the grid was found to be sufficient for the mean flow field to reach 
axisymmetric steady states at the computational inflow and outflow boundaries. There are 
81 axial points on the upper and lower blade surfaces, and 30 axial points on the upstream 
and downstream periodic boundaries. This distribution was found to be sufficient for most 
of the calculations reported herein, with approximately 20 points per wave being applied to 
resolve the dominant acoustic waves. However, for some of the 3D Tenth Standard Cascade 
calculations, the near-sonic conditions on blade suction surfaces resulted in short wavelength 
acoustic response phenomena that could not be resolved on the prescribed 141 x 41 x 11 
H-mesh. 

The TURBO near-field, finite-volume solutions have been coupled to far-field acoustic 
eigensolutions, which have been determined on a radial grid consisting of 24 points clustered 
near the hub and duct walls. For the present calculations, any nearly convected distur- 
bances that occur downstream of the blade row are simply convected numerically through 
the computational outflow boundary and into the fax downstream region of the flow. 

The full potential steady and the LINFLO linearized unsteady solutions were determined 
on composite meshes consisting of local C-meshes embedded in global H-meshes, which 
extended one axial chord upstream and downstream from the blade row. The H- and C- 
meshes used with LINFLO consisted of 155 axial and 41 tangential lines and 101 radial and 
21 circumferential lines, respectively. Coarser H- and C- meshes were used for the CASPOF 
calculations. 

The n um erical solutions, reported herein, were determined on an IBM-3CT Workstation. 
TURBO, “time-accurate,” steady, subsonic, inviscid solutions required 780 CPU minutes 
per 1,000 time steps and a minimum of 1,500 to 2,000 time steps to converge. The TURBO 
unsteady calculations were started from the steady solution, and performed using 500 time- 
steps per cycle of blade motion and four Newton iterations per time step. For single-passage 
solutions, six to eight cycles of motion were needed to converge the nonlinear inviscid so- 
lutions to a periodic state. The subsonic inviscid calculations required 350 CPU minutes 
per blade passage per cycle of blade motion. The number of blade passages included in a 
nonlinear unsteady calculation depends upon the interblade phase angle. For example, if 
a = 60 deg, six passages are needed. 

TURBO nonlinear steady and unsteady calculations require approximately 1350 /isec./ time- 
step /cell with four sub-iterations. The memory requirement, using 32-bit arithmetic, is ap- 
proximately 1.8 kilobytes /cell. This requirement is based on the option of using two blocks 
per blade passage and in-core storage for all variables. 

Steady Flow 

Predicted distributions of relative, steady, isentropic, surface, Mach number based on 
local static pressure [P(r, #,£)] and the local inlet relative total pressure [Pt,-oo(**)]> ie > 
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for the 3D and 2D, 10th Standard Cascades, axe shown in Figure 4. The inlet and exit, mean- 
flow quantities for the 3D calculation are given in Figure 5. For the CASPOF, full potential 
calculation, the relative inlet Mach number, M_ «, = 0.7, and inlet flow angle, ft-oo = 55 deg, 
are prescribed and a Kutta condition is imposed at the blade trailing edges. For the TURBO 
calculation, the total temperature (T| bs = 5.766) and total pressure (P£ bs = 2.237) of the 
uniform mean inlet flow (V 0 abs ) = 0 are specified at the computational inlet boundary (i.e., 
at f = — Cax), and the mean-flow static pressure at the hub is specified at the 

computational exit boundary (£ = 2cax), so that the relative inlet flow at midspan matches 
the 2D conditions. 

The TURBO steady-flow predictions at the hub, r/r D = 0.8, midspan, r/r D = 0.9, and 
tip, r/r D = 1.0, given in Figure 4, indicate that the Mach numbers on the blade suction and 
pressure surfaces show moderate variations with radius. Also, the 3D TURBO predictions 
at midspan are in close agreement with the 2D CASPOF predictions. The TURBO results 
indicate that the maximum Mach numbers on the suction surface of a blade are 0.849 at 
the hub, 0.906 at midspan, and 0.961 at the tip. These values occur at C/cax = 0.053, 0.073 
and 0.085, respectively. Thus, the flow is very close to sonic in the tip region, along a blade 
suction surface just aft of the leading edge. The CASPOF predictions for the 2D cascade 
indicate a maximum Mach number of 0.916 at £ = 0.065. 

For the three-dimensional flow, the steady static pressure (P = 1.4577), density ( p = 1.0), 
and axial velocity (V f = 0.5736) have constant values at inlet and the relative circumferential 
velocity, Vq = —Qr varies linearly from 0.7283 at the hub to 0.9103 at the tip. At the 
computational exit boundary, the steady pressure, density, and axial velocity vary with 
radius (mean shear), and the circumferential velocity varies nonlinearly with radius (mean 
swirl). As indicated in Figure 5, the steady blade loading causes increases in the pressure 
and density and decreases in the axial and circumferential velocities, especially the latter. 

Blade Vibration 

The 3D TURBO and LINFLUX analyses and the 2D LINFLO analysis have been applied 
to predict the unsteady aerodynamic responses of the 3D and 2D 10th Standard Cascades to 
pure bending and pure torsional blade vibrations at unit frequency, as described below. In 
a linearized analysis, such as LINFLUX or LINFLO, the far-field conditions are determined 
based on the input mean flow before the unsteady computation begins. In the nonlinear 
TURBO analysis, a converged solution, which is steady in the rotating frame of reference, is 
used as the initial solution for an unsteady computation. Temporal Fourier decompositions 
of the flow quantities at inlet and exit are performed as the solution is marched in time. The 
temporal Fourier coefficients are updated N times per cycle, where N is a user input. 

All of the results presented in this report have been determined using five updates per 
cycle. Thus, the far-field conditions are determined based on the most current temporal 
Fourier coefficients. Ideally, the zero-frequency, temporal Fourier coefficients should cor- 
respond to the initial mean flow. However, the back pressures used for all the unsteady 
TURBO calculations are the same as the back pressure used in the mean-flow calculation. 
As a result, the mean mass flow for different unsteady cases can differ by as much as 1% of 
the steady mass flow. One could maintain the same mean mass flow by adjusting the back 
pressure for each unsteady run, but this would be laborious and computationally expensive. 
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In the current 3D TURBO analysis, two to three temporal harmonics are kept in the 
far-field analysis. These higher harmonic terms, as will be shown later, are usually small 
compared to the first harmonic term. Thus, the classification of sub- or superresonant un- 
steady motion is based on the propagation properties of the fundamental (i.e., first-harmonic) 
acoustic disturbances. For example, the unsteady excitation at u = 1 and o = 90 deg is clas- 
sified as superresonant, because propagating acoustic response disturbances at the excitation 
frequency exist in the upstream region. If this excitation is a prescribed blade vibration, only 
acoustic response disturbances will occur in the far-field. 

Local ( wc ) and global (W c ) work-per-cycle predictions for the 3D 10th Standard Cascade 
undergoing pure torsional and pure bending vibrations at u = 1 and a = ±90 deg (No = ±6) 
are shown in Figure 6, where the TURBO local response predictions are given at eleven span- 
wise stations from hub (r/r D = 0.8) to tip (r/r D = 1.0). These results indicate that the 
local work per cycle responses to the blade torsional and bending blade vibrations do not 
vary significantly with radius, but, the results for the bending vibrations show greater radial 
variations than those for the torsional motions. Note that the multi-passage TURBO solu- 
tions show slight blade-to-blade variations. Hence, local work-per-cycle predictions shown 
in Figure 6 are those on the reference blade; and the global work-per-cycle predictions are 
the averaged values, taken over all blades operating within the numerical solution domain. 
We will discuss the blade-to-blade variation later in this section. 

The averaged local work per cycle predictions at midspan, as determined from the 3D 
TURBO, LINFLUX and the 2D LINFLO predictions, for the 10th Standard Cascade vi- 
brating in torsion and bending are shown in Figures 7 and 9, respectively, for blade mo- 
tions at unit frequency and at interblade phase angles, a, of —90 deg, 0 deg, +90 deg and 
+180 deg. The motions at a = Odeg and 90 deg are superresonant. For the in-phase mo- 
tions at a = Odeg, propagating acoustic response disturbances, at (m, //) = (0,0), occur 
both upstream and downstream of the blade row. For the motions at a = 90 deg, such a 
disturbance occurs only in the upstream region. For the (subresonant) motions at -90 deg 
and a = 180 deg all acoustic response disturbances attenuate. 

Figure 7 shows that the torsional response predictions determined using the TURBO code 
are in good agreement with the corresponding 3D LINFLUX and 2D LINFLO predictions. 
The TURBO calculations were run with radial-mode far-field boundary conditions. To 
demonstrate the need for such conditions, the same torsional vibration cases were run, with 
TURBO, using quasi-unsteady, local, one-dimensional, far-field boundary conditions [Jan89]. 
Figure 8 shows that the predictions based on the ID far-field conditions do not agree very well 
with the 2D LINFLO results, except for the subresonant vibration at a = -90 deg. When the 
blades undergo a subresonant motion, all acoustic response disturbances attenuate. Thus, 
ID boundary conditions are often adequate. However, for the vibration at a = 180 deg, 
which is also subresonant, the local work per cycle predictions show appreciable differences 
on the suction surface. Similar discrepancies are also noted for the a = 0 deg case. Although 
the unsteady motion at a = 0 deg is superresonant, ID boundary conditions should be 
capable of handling the planar propagating waves. Of the four cases shown in Figure 8, the 
unsteady predictions for the superresonant vibration at c = 90 deg clearly show the poorest 
agreement. 

Predictions for bending vibrations, as determined using the three aforementioned codes, 
are shown in Figure 9. Those for the bending vibrations at a = 0 and 180 deg show small 
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differences over the entire blade. The reasons for these differences are not understood at 
present, but similar discrepancies have been reported in earlier work in which the predictions 
of 2D nonlinear [AV94, AV96], 2D linearized (LINFLUX) [VMK95, MV95], and 3D linearized 
(LINFLUX) [MV97] analyses were compared with LINFLO results. Similar to the results 
of the 3D LINFLUX analysis, the local work per cycle predictions for the bending vibration 
at a = 90 deg show small differences along the pressure surface, but large differences on the 
suction surface. The reasons for the large discrepancies have not been established at this time. 
However, we suspect that local, high-wave- number, acoustic responses, occurring in regions 
of high-subsonic steady Mach number, are not adequately resolved on the 141 x 41 x 11 H- 
mesh used for the TURBO calculations. The local work per cycle predictions for the bending 
vibration at a = —90 deg also show important differences, in this case, on both the suction 
and pressure surfaces. This is the only TURBO run that is significantly different from the 
3D LINFLUX results given in [MV97]. Again, it is not understood why the discrepancies 
occur. 

Global work-per-cycle predictions for the 2D and 3D 10th Standard Cascade cascades 
undergoing prescribed blade vibrations are shown in Figure 10, where results for the global 
work per cycle versus interblade phase angle are given for pure torsional vibrations about 
midchord and pure bending vibrations at unit frequency. The 3D TURBO results, indicated 
by the circular symbols in Figure 10, have been determined for N D = -6, -4, 0, 4, 6, 8, 16, 
and 18; the 3D LINFLUX results, by the square symbols, for N D = -6, -5, . . . , 18, and the 
2D LINFLO results, for -90 deg < cr < 270 deg in increments of one degree. The 2D work 
per cycle predictions are multiplied by the blade span, i.e., td - r H = 0.2 r D = 0.849, to 
allow a convenient comparison with the 3D predictions. 

The resonance or cut-off conditions for the two-dimensional configuration are aZ^ = 
-26.93 deg and = 117.12 deg in the far upstream region and = -31.80 deg and 
< 7 +^ = 59.79 deg in the far downstream region. The superresonant blade motions at u = 
1 occur at interblade phase angles between these cut-off values and send a propagating 
wave into the upstream and/or downstream regions of the flow. The blade motions at 
-90 deg < a < -31.80 deg and 117.12 deg < a < 270 deg are subresonant. The results in 
Figure 10 indicate a very good agreement between the 3D TURBO, the 3D LINFLUX and 
the 2D LINFLO global response predictions over the entire nodal diameter or interblade 
phase angle range of blade vibrations. We should reiterate, however, that for superresonant 
bending vibrations at a = 90 deg, and the subresonant bending motions at a = -90 deg, 
the TURBO and LINFLO local responses show large differences, cf. Figure 9. 

Next, we consider the acoustic properties far from the blade row. Predicted steady, 
as well as first- and second-harmonic axial eigenvalues and first-harmonic radial pressure 
modes, p^(r), for m = -1,0,1 and p = 0,1 axe shown in Figures 11 through 13. Here, 
the unsteady excitation occurs at u = 1 and Nd = 6 (a = 90 deg). Because of mean blade 
l oading , the steady inlet and exit conditions for the 3D 10th Standard Cascade differ. As 
a result, the acoustic properties, Xm» and in the far-upstream region of the flow, differ 
from those in the far-downstream region. In particular, for an unsteady excitation at u = 1 
and a = 90 deg, the fundamental acoustic disturbances in the (0,0) mode are of propagating 
type fax upstream, but, of attenuating type far downstream. Note that the second-harmonic 
acoustic disturbances in the (0,0) mode axe of propagating type in both the upstream and 
downstream regions. 
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Tn the far-upstream region of an unsteady flow at lo — 1 and o — 90 deg through the 3D 
10th Standard Cascade , the propagating acoustic response disturbance has an axial wave 
number, of 1.584 and the least damped or (0,1) response disturbance has an attenuation 
constant, /?, of 3.964. In the far-downstream region, p = -1.077 for the least-damped or 
(0,0) response disturbance. These numbers agree favorably with those predicted by the 
3D LINFLUX run, namely, 1.583, 3.990, and -1.084. Since the absolute far-downstream 
mean flow is nonuniform, the axial wave numbers of the attenuating disturbances in a given 
circumferential mode vary with radial mode number, p, as indicated in Figure 12, particularly 
those for m = 1. 

The radial eigenmodes for the pressures associated with the far upstream acoustic exci- 
tations or responses and the far downstream acoustic responses for an unsteady excitation 
at u = 1 and o — 90 deg are shown in Figure 13. Although the inlet and exit mean-flow 
conditions differ, the upstream and downstream radial pressure modes are very similar, with 
the downstream modes showing a somewhat greater radial variations than their upstream 
counterparts. Note that the phase of a modal pressure disturbance is independent of ra- 
dius for the uniform absolute mean flow at inlet, but the phase varies with radius for the 
mean flow with swirl and axial shear that exists in the far downstream region. Thus, the 
far-upstream, pressure modes, p R (r), are purely real, but the far-downstream modes have 
some imagin ary or out-of phase content. 

The TURBO calculations for the subresonant a = -90 deg and the superresonant a = 
90 deg blade motions reveal that, for the most part, the far-field acoustic responses are of 
small amplitude at the computational inflow and outflow boundaries. However, the super- 
resonant torsional and bending vibrations at o = 90 deg produce upstream propagating 
acoustic response disturbances which have amplitudes, a a, of 1.352 and 1.540, respectively, 
and occur at an axial wave number, of 1.583. The corresponding LINFLO predictions 
are a A = 1.529 and 2.822 and = 1.603. Thus, there is a substantial difference between 
the TURBO and LINFLO predictions for the upstream propagating, (0,0), acoustic response 
waves caused by the bending vibration. 

In Table 6.1, the TURBO-predicted eigenvalues and amplitudes of the upstream prop- 
agating, (0,0), acoustic response wave caused by the bending and torsion vibrations are 
compared ag ains t the corresponding LINFLUX and LINFLO values. Note that TURBO 
updates the far-field eigenmodes several times per cycle. Certain eigenmodes, usually higher 
order ones, may be missing from one update but reappear in the next update because the 
current eigenmode filtering scheme sometimes filters out modes that should be retained. In 
Table 6.1 the eigenvalue for the bending vibration is slightly different from the one for the 
torsional vibration, because the meanflow evolves differently in the two unsteady compu- 
tations. Also, the eigenvalues and amplitudes for the (0,0) acoustic response mode of the 
second temporal harmonic are listed in the table. The second-harmonic content is significant 
in the initial cycle, but becomes negligible after eight cycles. 

Next, we examine the results for a bending vibration at u = 1 and a = —90 deg in some 
detail. This case is chosen because the nonlinear response predictions are in poor agreement 
with the results of both 2D and 3D linearized analyses. In addition, the response predictions 
show the most blade-to-blade variation. Furthermore, this is the only case that we were 
not able to get a converged work-per-cycle prediction even after more than forty cycles of 
unsteady computations. 


30 





dA 

torsion 

bending 

torsion 

bending 

LINFLO 

1.603 

1.529 

2.822 

LINFLUX 

1.583 

1.352 

1.540 

TURBO (1st harmonic, 8th cycle) 

1.587 

1.584 

1.375 

1.610 ] 

TURBO (2nd harmonic, 1st cycle) 

3.216 

3.218 

0.103 

0.103 

TURBO (2nd harmonic, 8th cycle) 

3.223 

3.219 

5.99 x 10" 3 

4.96 x 10" 4 


Table 6.1: Comparison between the axial wave numbers (k^) and acoustic disturbance ampli- 
tudes (oa) predicted by the 2D linear analysis (LINFLO), the 3D linear analysis (LINFLUX), 
and 3D nonlinear analysis (TURBO). 

In Figure 14, the calculated time histories of the mass flow, m = f pV • dA, and total 
pressure ratio, Pt,+ 00 / Ft ,— 00 for the first twenty cycles are shown. It can be observed that 
the mass flow reached a periodic state in approximately six cycles. The normalized mean 
mass flow of 8.125 is approximately 0.7% lower than the normalized steady mass flow of 8.18. 
The time history of the mass flow shows significant higher harmonic content, even after the 
mass flow reaches a periodic state. We have not been able to determine the cause of the high 
non-harmonic variation in mass flow that occurs during the initial cycles. The inlet/exit mass 
flow fluctuation is observed immediately after the start of the unsteady computation, thus, 
this fluctuation must be generated in the far field, not in the interior of the solution domain 
due to the blade motion. One possible source of numerical error might be a deforming grid 
formulation, in which the “geometric conservation law” is violated. However, to the best 
of our knowledge, the TURBO deforming grid formulation and implementation are correct 
and do not violate the geometric conservation law. Moreover, a similar implementation was 
successfully used in the 2D nonlinear NPHASE analysis [AV94, AV96]. 

In the current TURBO unsteady analysis, we have included two to three temporal har- 
monic terms in the far-field analysis to minimize the reflection of outgoing transient higher 
harmonic waves back into the solution domain. Unfortunately, in addition to the afore- 
mentioned non-harmonic terms appearing in the far-field, the initial impulsive blade motion 
appears to have generated anharmonic and higher harmonic waves. Some of these waves are 
reflected, and the reflections ultimately result in blade-to-blade variations in the work-per- 
cycle response. This effect seems to be most noticeable for cases in which blades undergo 
bending vibrations. 

Detailed response results for the 3D 10th Standard Cascade undergoing a bending vi- 
bration at u = 1 and a — -90 deg are shown in Figures 15-17. An examination of the 
time-mean static pressure shows that the time-mean loadings on all four blades are the 
same. The real and imaginary parts of the first-harmonic unsteady pressure at mid-span 
are shown in Figure 15. The unsteady pressures on the neighboring blades do appear to be 
90 deg out of phase. For example, the in-phase unsteady pressure on blade 0 (solid line) is 
similar to the out-of-phase unsteady pressure on blade 1 (long dashes). However, the local 
work-per-cycle predictions in Figure 16 show noticeable variations from blade to blade. The 
midpsan, local work-per-cycle predictions on all the four blades are shown in Figure 17. For 
comparison with the results of the 2D linearized analysis, we have averaged the TURBO, 
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midspan, local work-per-cycle predictions over all blades, cf. Figures 7 and 9. 

6.2 NASA Rotor 67 

The NASA Rotor 67 is a research transonic fan consisting of 22 blades. The tip diameter 
of the fan rotor varies from 51.4 cm at the leading edge to 48.5 cm at the trailing edge, and 
the hub-to-tip radius ratio varies from 0.375 to 0.478 from the inlet to exit. At the design 
point, the rotational speed of the rotor is at 16,043 rpm. With an inlet axial Mach number 
of approximately 0.49, the tip speed at the design condition is 429 m/sec, which corresponds 
to a tip relative Mach number of 1.38. Also, the mass flow rate at the design point is 33.25 
kg/sec and the total pressure ratio is 1 .63. 

The purpose of the present numerical study is to demonstrate the capabilities of the 
TURBO analysis for analyzing the flutter characteristics of realistic transonic fans. Our 
goal is to obtain steady inviscid solutions at two throttle positions on the design speed line, 
one near peak efficiency and the other near stall. At each point, we will perform flutter 
analyses for a bending and a torsional vibration at two different nodal diameters. 

The grid for the Rotor 67 calculations, see Figure 18, consists of 121 axial, 33 tangential 
and 17 radial surfaces (61,440 cells), and extends one axial chord, at mid-span, upstream and 
downstream from the blade row. For the calculations, we assume that there is no clearance 
between the rotor blades and the outer duct wall. We have also idealized the endwalls such 
that near the computational inlet and exit boundaries, the inner and outer duct radii are 
constant. This is necessary because the 3D far-field eigensolver assumes that the inlet and 
exit flows are fully developed and, therefore, that they do not vary with axial distance. The 
computational grid consists of 65 axial points on the upper and lower blade surfaces, and 
28 axial points on the upstream and downstream periodic boundaries. The grid is clustered 
near blade leading and trailing edges, and near blade surfaces. The current version of the 
TIGER grid generation package [SS91] used for generating the Rotor 67 grid, does not have 
an elliptic grid smoother. Thus, the grid quality near the leading- and trailing-edge planes is 
relatively poor. Furthermore, we have not examined whether the selected grid has adequate 
resolution for the Rotor 67 flows being considered, especially in the axial direction. The size 
of the grid was chosen so that the present unsteady calculations could be performed on a 
128MB workstation. 

TURBO time-accurate steady subsonic inviscid solutions for the Rotor 67 study required 
680 CPU minutes per 1,000 time steps on an IBM-3CT Workstation. It is unclear to us why 
this case with slightly more grid cells (61,440) than that of the 3D 10th Standard Cascade case 
(56,000) require less CPU time per 1,000 time steps. However, for the Rotor 67 analysis, 
the first converged steady-state solution took more thai 8,000 time steps. The TURBO 
unsteady calculations were started from the appropriate steady solution, and performed 
using 500 time-steps per cycle of blade motion, four Newton iterations per time step. For 
each cycle of blade motion, the unsteady inviscid calculat ions required 375 CPU minutes per 
blade passage. 

As for the calculations for the 3D 10th Standard Cascade, the TURBO near-field, finite- 
volume solutions are coupled to far-field acoustic eigensolutions, which are determined on 
a radial grid consisting of 24 points clustered near the hub and duct walls. Any nearly 
convected disturbances that occur downstream of the blade row are simply convected nu- 
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merically through the computational outflow boundary. 

Steady Flow 

Numerous researchers; e.g., Chima [Chi91], Jennions, et al. [JT92], Amone [Arn93], Rhie, 
et al. [RZH+93], have performed numerical simulations of viscous flows through the NASA 
Rotor 67. Most showed the flow field to be very complicated and characterized by effects 
such as shock-boundary layer interaction, tip-leakage flow, unsteady vortex shedding, and 
flow separation with vortex roll up. Under the present effort, we have performed mviscid 
simulations of steady and unsteady flows. 

Initially, we experienced difficulties in getting converged steady-state inviscid solutions. 
Consequently, we had to use very small time steps for the first quarter of the rotor rotation 
before the time step or the CFL number could be increased. Also, because of the relatively 
coarse grid used for our calculations, we did not expect the resulting solutions to agree well 
with data. As a result, we did not run the whole speed line. Instead, we ran two different back 
pressures to steady state to approximate the peak efficiency and near stall points reported 
in [SWHS89]. The mass flow for the first point (near peak efficiency) is 34.8 kg/sec and the 
total pressure ratio is 1.67. The mass flow for the second point (near stall) is 31.0 kg/sec 
and the total pressure ratio is 1.75. As expected, the inviscid “speedline is higher than the 
experimental speedline due to lower losses, see Figure 19. Surprisingly, the TURBO inviscid 
solutions show good qualitative agreement with the experimental data. 

In Figures 20 and 21, the numerical and experimental relative Mach number contours 
at three different spanwise locations are shown for the two different operating conditions 
mentioned above. Note, the constant J surfaces on which the contours plots of the numerical 
solution are shown are not constant-radius surfaces. Also, the contour plots of the numerical 
solutions are plotted in the f , 0-plane instead of the £,r 0-plane because the plotting package 
used in generating the contour curves could not properly handle periodic boundaries that 
have varying radial locations. This explains why the airfoil geometries in the contours plots 
for the experimental data are not the same as those in the contour plots for the numerical 

results. . 

The Mach number contours near the peak efficiency point are shown in Figure 20. Those 

at 10 percent of span from the tip show an inlet Mach number of 1.35. There is a bow shock 
at the leading edge, and a normal in-passage shock near the trailing edge on the suction 
surface. At 30 percent of span, the flow pattern is similar, except that the inlet Mach 
number is lower, i.e., approximately 1.2. At 70 percent of span, the inlet Mach number is 
0.95 and a supersonic bubble appears on the suction surface near the leading edge. 

In Figure 21, we consider the near stall condition. The Mach contours at 10 and 30 
percent of span ’show that the location of the in-passage shock is near midchord on the 
suction surface, and the strength of the shock is much stronger than the normal shock that 
occurs at near peak efficiency. Thus, the exit Mach number is lower in the near stall case. 
At 70 percent span, the size of the supersonic bubble has increased significantly. In general, 
the qualitative agreement between experiment and the computations for both operating 
conditions is good. However, because of grid clustering near blade edges and near blade 
surfaces, the grid is relatively coarse in the mid-chord and mid-passage regions. Hence, the 
inviscid solutions show highly smeared shocks. 
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Figure 22 shows a carpet plot of the chordwise steady pressure distributions at 17 span- 
wise stations. Here, the axial distance is normalized by the axial chord length at the hub. 
Due to blade twist, the normalized axial locations of the leading and trailing edges at the 
blade tip are approximately 0.25 and 0.75, respectively. The curves in Figure 22 clearly 
depict the differences between in-passage shock locations for the two steady solutions. These 
curves have pressure spikes at the blade trailing edge that are due to flow overspeeds over 
blunt trailing edges. 

The radial distributions of the inlet and exit mean-flow quantities for Rotor 67 operating 
near peak efficiency are shown in Figure 23. At inlet, the relative circumferential velocity, 
V e = -Dr varies linearly from 0.436 at the hub to 1.263 at the tip, and the steady static 
pressure (P), density (p), and axial velocity (V f ) are nearly constant. At the computational 
exit boundary, the steady pressure, density, and axial velocity vary with radius (mean shear), 
and the circumferential velocity varies nonlinearly with radius (mean swirl). As indicated in 
Figure 23, the steady blade loading causes increases in the pressure and density and decreases 
in the circumferential velocity. There is only a slight decrease in the axial velocity because 
the steady-flow operating point is not too far from the choke point. 

Blade Vibration 

For the flutter analysis, we have arbitrarily chosen a blade vibratory frequency of 1.19 
times the rotation speed, giving a reduced frequency of 0.54 (based on the midspan relative 
inlet velocity and blade chord), because we do not have any structural information for the 
blade. This frequency is probably representative of, or slightly higher than, the first bending 
mode frequency of a typical low aspect ratio fan. Presently, a torsional and bending mode of 
the same frequency have been studied in our flutter analysis of Rotor 67. As stated earlier, 
simple 2D analytical mode shapes have been used, because there is no readily available 
structural information for the Rotor 67 configuration. We should note that the TURBO code 
is capable of handling finite-element mode shapes, such as NASTRAN-generated modes. 

The TURBO analysis was run for vibrations at two nodal diameters, N D = 0 and N D = 
11, for each of the meanflow conditions discussed above. At the near peak efficiency point, 
the predicted mean-flow, first-harmonic, and second-harmonic axial eigenvalues, at the inlet 
and exit, are shown in Figures 24 and 25, respectively, for the acoustic modes at m = 
-1,0,1, n = 0,1, 2, 3, 4 and an unsteady excitation at u = 0.54 and N D = 0 (a = Odeg). 
The first-harmonic radial pressure modes, Pmft( r )> 8X6 shown in Figure 26. As for the subsonic 
3D 10th Standard Cascade, the steady inlet and exit conditions differ due to mean blade 
loading. Thus, the acoustic properties and p£ M , in the far-upstream region of the flow 
differ from those in the far-downstream region. However, unlike the results for the subsonic 
3D 10th Standard cascade, there are many propagating sicoustic modes in the far-upstream 
region because the meanflow at inlet is supersonic in the blade-tip region. 

In the far-upstream region of an unsteady flow, at i v = 0.54 and a = Odeg, through 
Rotor 67 operating near peak efficiency, the far-field etgenanalyses of both the first and 
second temporal harmonics show that the first three radial modes for m — 1, the first two 

modes for m = 1, and the (0,0) mode are all propagating. For the first temporal harmonic, 
the axial wave number, k^, for modes (-1,0), (0,0), and (1,0) are 111.8, 6.766, and -94.4, 
respectively. As a result, the computational grid used in the near-field calculation does not 
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have the resolution needed to accurately resolve the (-1,0) and (1,0) modes. For the torsional 
vibration, The least-damped or (0,1) response disturbance has an attenuation constant, /?, 
of 10.86. The (0,0) upstream propagating acoustic response disturbance has an amplitude, 
a A , of 1.35 which is linearly scaled to correspond a blade pitching amplitude of one radian 
about the midchord locations. All the other modes are at least two orders of magnitude 
smaller than the (0,0) acoustic response disturbance. However, near stall, the amplitudes 
of the (-1,2) and (1,1) acoustic disturbances are approximately 20% and 10% of that of the 
(0,0) acoustic disturbance. 

Far fewer modes are determined in the far-downstream region than in the far-upstream 
region. This is probably an artifact of the numerical filtering scheme currently used in the 
far-field eigenanalysis. For this particular case, the filtering scheme did not eliminate all the 
superfluous modes. For example, there is one superfluous mode for the first temporal har- 
monic. We do not know whether these nonphysical modes will degrade the solutions or not. 
Only the (0,0) acoustic response disturbances are propagating. For the torsional motion, the 
first-harmonic axial wave number of the (0,0) mode is 5.88, and the amplitude of the distur- 
bance is 0.123. The least-damped or (0,1) response disturbance has an attenuation constant 
of -15.82, and a non-negligible amplitude of 0.058. For the second temporal harmonic, the 
axial wave number of the (0,0) mode is 11.73, and the amplitude of the disturbance is 0.103. 
Also, the attenuation constant of the (0,1) mode is -14.48, and the amplitude of this mode 
is 0.032. 

The radial eigenmodes for the pressures associated with the far upstream acoustic exci- 
tations or responses and the far downstream acoustic responses for an unsteady excitation 
at co = 0.54 and a = 0 deg are shown in Figure 26. In addition to the differences between the 
inlet and exit mean-flow fluid properties, the hub-to-tip ratios are different at inlet and exit. 
Yet, the lowest order upstream and downstream radial pressure modes are very similar. 

Local (w c ) and global (W c ) work per cycle predictions for Rotor 67 undergoing pure 
torsional motions at u = 0.54, and <7 = 0 and 180 deg {Nd = 0 and 11) at the two aforemen- 
tioned meanflow conditions are shown in Figure 27. The TURBO analysis indicates that the 
blade motions are stable. Moreover, as indicated by the global work per cycle predictions, 
the blade motions at N& = 11 axe more stable than those at N D = 0. Also, the results 
show the blade to be more stable when operating at the near stall point than at the near 
peak efficiency point. The last observation seems to contradict previous experience. A close 
examination of the local work per cycle distributions shows the unsteady forces around the 
in-passage shock play an important role in determining stability of the blade motion. Near 
the peak efficiency point, the location of the in-passage shock in the outer span region of the 
blade is near the trailing edge on the suction surface and near mid-chord on the pressure 
surface. At N D = 0, the unsteady forces ahead of the shock on both the pressure and suc- 
tion surfaces extract energy from the torsional blade motion. For an out-of-phase torsional 
motion, the unsteady forces on the pressure surface ahead of the shock are significantly more 
stabilizing than are those for the N D = 0 case, whereas forces ahead of the shock on the 
suction surface are destabilizing. Near stall, the in-passage shock moves forward towards 
the leading edge on the pressure surface, and towards mid-chord on the suction surface. 
The unsteady forces on the pressure surface near the leading edge, where the shock sits, are 
extracting work from the blade; whereas the unsteady forces on the suction surface in the 
vicinity of the shock, are stabilizing for N D = 0 and destabilizing for N D = 11. 
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Local and global work per cycle predictions for Rotor 67 undergoing pure bending motions 
at u = 0.54, o = 0 and 180 deg and for the two meanflcw conditions being considered are 
shown in Figure 28. Once again, the TURJBO predictions show the blade motion to be stable 
at the four points analyzed. As for the torsional vibrations, the blade motions are more stable 
at N d = 11 than at Nd = 0. However, as the back pressure increases, the bending vibrations 
at cr = 0 deg become less stable. Again, as for the torsional vibrations, the unsteady forces 
near the shock dominate the work-per-cycle predictions for the peak efficiency point. At 
this meanflow condition, the unsteady forces on the pressure surface ahead of the shock axe 
stabilizing, while the unsteady forces on the suction surface are destabilizing, if N D = 0 and 
stabilizing, if N D = 11. At the near stall meanflow condition, the contributions to the local 
work per cycle from the unsteady loads away from the shock are significant, especially for 
the bending vibration at N D = 0, in which the unsteady forces over the latter half of the 
blade on the suction surface, contribute energy to the blade motion. 

6.3 Discussion 

We have presented numerical results for unsteady flows through a three-dimensional ver- 
sion of the 10th Standard Cascade. These results pertain to flows in which the unsteady 
fluctuations are caused by prescribed blade vibrations. They were determined using the 3D 
TURBO and LINFLUX analyses and the 2D LINFLO analysis. TURBO employs an im- 
plicit, flux-split, finite-volume scheme for solving the unsteady Euler equations in the near 
field, which typically extends from one axial chord upstream to one axial chord downstream 
of the blade row, and numerical eigenanalyses for determining unsteady perturbations of 
fully-developed, axisymmetric, swirling mean flows in the far upstream and far downstream 
regions. For the unsteady flows considered herein, the eigenanalyses have been used to 
determine the first two or three modal acoustic disturbances, additional higher-order distur- 
bances are ass um ed to be of negligible amplitude at the computational inflow and outflow 
boundaries, and the remaining part of the unsteady perturbation, consisting of convected 
and nearly convected disturbances, is simply convected out of the near-field domain through 
the computational outflow boundary. 

The numerical results indicate that the far-field eigenanalysis is capable of providing 
reasonable solutions for the axial eigenvalues and the radial pressure modes (e.g.,see Fig- 
ures 11, 12, and 13) of the acoustic excitations and responses that can exist far upstream 
and far downstream of a blade row. At this point, we have not applied the eigenanalysis 
to predict the axial eigenvalues and radial eigenmodes associated with nearly convected dis- 
turbances. The behavior of such disturbances is not well understood at present, as far-field 
eigenanalyses for non-uniform mean flows have become available only recently. However, 
it will be necessary to provide accurate numerical representations of nearly-convected, pre- 
dominantly vortical, disturbances to predict the unsteady aerodynamic responses associated 
with wake/blade-row interactions. 

The TURBO predictions for the zeroth-order or steady relative flow at _ <*, = 

0.4015 through the 3D 10th Standard Cascade shows moderate variations in the blade-surface 
Mach numbers with radius, (Figure 4), and small variations in blade loading. In addition, the 
3D Euler predictions for the surface Mach numbers at blade midspan are in close agreement 
with 2D full-potential predictions. The 3D 10th Standard Cascade operates in a uniform, 
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axial, absolute, mean inlet flow, but, because of steady blade loading, the mean flow far 
downstream of the blade row (Figure 5) has swirl and axial shear. 

The TURBO local unsteady response predictions , i.e., wc vs f, (Figures 6, 7 and 9) for 
the 3D 10th Standard Cascade, undergoing pure torsional and pure bending vibrations at 
u = 1, show small variations with radius and, for the most part, the results at midspan, 
are in good agreement with the 3D LINFLUX and 2D LINFLO predictions. However, the 
TURBO and LINFLUX local work-per-cycle results for a superresonant bending vibration at 
a = 90 deg are significantly different from the 2D LINFLO predictions. We suspect that these 
differences are due to an inadequate resolution, of the local, high wave number, upstream 
traveling, acoustic response disturbances that occur at high-subsonic Mach numbers. For the 
subresonant bending vibration at o — —90 deg, the 3D LINFLUX and 2D LINFLO results 
show good agreement, but the TURBO predictions do not agree with those of the two linear 
analyses. Furthermore, the convergence of the work-per-cycle for this subresonant bending 
vibration case using TURBO is very slow. We do not have an explanation for the slow 
convergence of the work-per-cycle calculations, and the discrepancies between the TURBO 
and LINFLUX solutions along the blade surfaces. 

The TURBO, LINFLUX and LINFLO global work per cycle, W c vs a, predictions for 
torsional and bending vibrations (Figure 10) are in very good agreement. However, the 
global results for the bending vibrations must be interpreted with some caution, as the local 
responses differ along surfaces at several interblade phase angles. 

We have presented some numerical results for unsteady flows through a realistic tran- 
sonic fan, i.e., the NASA Rotor 67. Our ultimate goal is to use the 3D TURBO analysis 
for aeromechanical stability assessment of turbomachinery blades. For the unsteady flows 
considered herein, the eigenanalyses have been used to determine the first two or three modal 
acoustic disturbances. Although, we have a reasonable understanding of the eigenanlysis for 
subsonic mean flows, such as those associated with the 3D 10th Standard Cascade, we do 
not know how reliable this analysis is for supersonic flows. As discussed previously, some- 
times, there axe missing or superfluous modes. We do not know whether or not these modes 
seriously deteriorate the solutions. 

Nevertheless, some interesting trends have been observed from the results of the TURBO 
unsteady analysis of Rotor 67 undergoing bending and torsional vibrations. The analysis 
shows the bending vibrations at Nd — 0 and Nd = 11 to be quite stable, much more so than 
for the blades undergoing corresponding torsional vibrations. Also, for both torsional and 
bending vibrations, the blades are less stable at Nd == 0 than at Nd = H- lu addition, the 
torsional mode of vibration becomes slightly more stable as back pressure increases, while 
the bending mode of vibration becomes less stable as back pressure increases. In general, 
all these trends seem to agree with those observed for classical supersonic unstalled flutter, 
which is usually associated with a torsional mode, and subsonic/transonic high-incidence 
flutter, which is usually associated with a bending mode. 
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7. Concluding Remarks 


The TURBO [JHW92], nonlinear, unsteady, aerodynamic analysis has been extended 
for turbomachinery aeroelastic applications. This analysis is based on the Euler/Navier- 
Stokes equations of fluid motion; a near-field, implicit, flux-split, finite-volume, analysis; 
and far-field eigenanalyses for the unsteady perturbations of fully-developed, axisymmetric, 
swirling mean flows. The far-field eigenanalyses, which are coupled to the near-field finite- 
volume analysis at computational inflow and outflow boundaries, allow incoming external 
aerodynamic excitations to be prescribed, and acoustic response disturbances to pass through 
these computational boundaries without spurious reflections. Under the current effort, no 
external aerodynamic excitations have been considered. 

We have applied the TURBO analysis to predict unsteady subsonic flows through a 
simple turbomachinery configuration, i.e., a three-dimensional version of the 10th Standard 
Cascade Configuration. We have also applied this analysis to predict unsteady flows through 
a realistic transonic fan, i.e., the NASA Rotor 67. We have considered unsteady flows excited 
by prescribed blade vibrations that are highly two dimensional. For the 3D 10th Standard 
Cascade, this allows us to compare and validate TURBO results against predictions based 
on previous two-dimensional analyses. 

The n um erical results indicate that the current version of the TURBO code is capable of 
providing accurate aerodynamic response information for unsteady subsonic flows, provided 
that the grids employed have a sufficient overall density and local clusterings in regions of 
high flow gradients. In particular, the numerical results indicate that the axial eigenvalues 
and radial eigenmodes of far-field acoustic disturbances can be accurately represented, and 
that the 3D blade-surface, response predictions show reasonable radial trends. The TURBO 
results at blade midspan and the 2D LINFLO results for the 3D 10th Standard Cascade axe 
in good qualitative agreement, but in some cases significant quantitative differences occur. 
The differences occur along the suction surfaces of the blades, where steady Mach numbers 
are close to one, and upstream of the blade row. Some evidence [MV97] suggests that 
the quantitative differences between the TURBO and LINFLO results can be eliminated if 
the meshes used in the TURBO calculations are of sufficient density and the grid lines are 
properly distributed. 

Based on the numerical results presented in this report, it appears that the far-field eigen- 
solver, developed for the TURBO code, is working properly and that it has been successfully 
coupled with the near-field numerical algorithm. Also, the TURBO analysis can yield useful 
response information for unsteady flows excited by blade vibrations. However, the proper 
treatment of superharmonic and anharmonic waves needs further investigation. In addition, 
the mesh requirements for accurately resolving such flows must be better understood. The 
requirements for flutter applications, for which reduced frequencies are typically of order 
one, can be readily met, but those for forced response studies, in which reduced frequencies 
on the order of 5 to 50 must be considered, will impose serious constraints on available 
computational resources. 

To improve the efficiency of TURBO steady and unsteady calculations, a number of 
computational strategies should be investigated. For example, a rotating-frame version of 
the TURBO analysis could be constructed to allow more efficient predictions of nonlinear, 
steady flows via the use of convergence accelerating schemes. Also, second-order accurate, 
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surface boundary conditions could be incorporated into TURBO to reduce the time required 
to achieve converged low-loss, steady and unsteady solutions. In addition, a single-passage 
version of the TURBO analysis with time-lagged periodic boundary conditions should be 
considered. This will greatly reduce the computational requirements for unsteady flows in 
which the structural vibratory pattern has a non-zero nodal diameter. Finally, a parallel ver- 
sion of TURBO should be considered, particularly for viscous and high-frequency unsteady 
flows. 

To date, we have focused on demonstrating and validating inviscid version of the TURBO 
code for flutter applications. In addtion, we have applied the code to predict flutter in 
transonic inviscid flows. However, the code needs further validation for transonic flows. Also, 
viscous effects are expected to play an important role in the analysis of high-incidence flutter. 
Thus, a validation the viscous capabilities of the TURBO code for unsteady applications 
should be carried out. 
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Figure 1: Rotating axial compressor blade row operating within an annular duct. 
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Figure 2: 3D Tenth Standard Configuration undergoing an exaggerated torsional motion 
(<*hub = 0 deg, a tip = 45 deg). The rotor consists of 24 airfoils. The nodal diameter of the 
blade motion is 6, which results in an interblade phase angle of 90 deg. The outer casing 
has been eliminated from the figure for clarity. 



Figure 3: TURBO computational grid at midspan for the 3D 10th Standard Cascade. 
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Figure 4: Relative steady isentropic surface mach number distributions for the 3D 10th 
Standard Cascade = 0.4015, |Q| = 0.2145). 
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10th Standard Cascade (Mi* = 0.4015, |Q| = 0.2145). 
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Figure 6: Local work per cycle distributions at all spanwise stations and global works per 
cycle for the 3D 10th Standard Cascade undergoing pure torsional vibrations about midchord 
and pure bending vibrations at u = 1 and u — ^90 deg (N& = ^6). 
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Figure 7: Averaged local work per cycle distributions at midspan, as predicted using the 
3D TURBO and the 2D LINFLO analyses, for the 3D 10th Standard Cascade undergoing 
torsional blade vibrations about midchord at a> = 1 . 
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Figure 8: Averaged local work per cycle distributions at midspan, as predicted using the 
3D TURBO analysis with one-dimensional, local, far-field conditions and the 2D LINFLO 
analysis, for the 3D 10th Standard Cascade undergoing torsional blade vibrations about 
midchord at u = 1. 
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Figure 9: Averaged local work per cycle distributions at midspan, as predicted using the 
3D TURBO and the 2D LINFLO analyses, for the 3D 10th Standard Cascade undergoing 
bending vibrations at u) = 1. 
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Figure 10: Work per cycle versus interblade phase angle for the 3D 10th Standard Cascade 
undergoing pure torsional vibrations about midchord (top) and pure bending vibrations 
(bottom) at u = 1. 
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Figure 11: Axial eigenvalues, x = P + at inlet for three circumferential (m = —1,0, 1) 
and three radial (/z = 0, 1,2) modes of steady (n = 0), iirst-harmonic (n = 1) and second- 
harmonic ( n = 2) acoustic disturbance in an unsteady flow at u = 1.0 and N D = 6, through 
the 3D 10th Standard Cascade. 
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Figure 12: Axial eigenvalues, x = P + ***> at exit for three circumferential (m = -1,0, 1) 
and three radial (fj, = 0,1,2) modes of steady (n = 0), first-harmonic ( n = 1) and second- 
harmonic (n = 2) acoustic disturbance in an unsteady flow at a; = 1.0 and N D = 6, through 
the 3D 10th Standard Cascade. 




Figure 13: First harmonic, radial, pressure modes, (r), m = —1,0,1, p = 0,1,2, at 
inlet and exit due to an acoustic excitation or response, at Nd = 6, far upstream and for 
an acoustic response, at N& = 6, far downstream of the 3D 10th Standard Cascade: ( — ) 
in-phase component; ( ) out-of-phase component of p^(r). 













Figure 15: First-harmonic unsteady pressure distributions at midspan acting on four blades 
of the 3D 10th Standard Cascade undergoing a pure bending vibration at a; = 1.0 and 
N d = -6. 
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Figure 16: Local work per cycle distributions on four blades of the 3D 10th Standard Cascade 
undergoing a pure bending vibration at u; = 1.0 and N D = —6. 
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Figure 17: Local work per cycle distributions at midspan for the 3D 10th Standard Cascade 
undergoing a pure bending vibration at u = 1.0 and Np = —6. 
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Figure 18: TURBO computational grid for the NASA Rotor 67 and leading edge detail at 
blade tip. 
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Figure 19: Comparison of measured and computed fan total pressure characteristics for the 
NASA Rotor 67 fan. 
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Figure 22: Steady static pressure distributions for the Rotor 67 fan (a) near the peak effi- 
ciency point, and (b) near stall. 
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Figure 24: Axial eigenvalues, x = P + at inlet for the time-mean (n = 0), and the first 
(n = 1) and second (n = 2) temporal harmonics, and three circumferential (m = —1,0, 1) 
and five radial (fj- = 0, 1, 2, 3, 4) modes of acoustic disturbance in the far field of Rotor 67, for 
an unsteady flow at u = 0.54 and N D = 0, and a meanflow condition near peak efficiency. 
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Figure 25: Axial eigenvalues, x = P + *«?> at exit for the time-mean (n = 0), and the first 
(n = 1) and second (n = 2) temporal harmonics, and three circumferential (m = -1, 0, 1) 
and five radial (j /u = 0, 1, 2, 3, 4) modes of acoustic disturbance in the far field of Rotor 67, for 
an unsteady flow at u = 0.54 and N D = 0, and a meanflow condition near peak efficiency. 
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Figure 26: First harmonic, radial pressure modes, p* M (7'), m = —1,0, 1, n = 0, 1,2, for an 
acoustic excitation or response, at No = 0, far upstream and for an acoustic response, at 
N d = 0, far downstream of the Rotor 67 fan: ( — ■) in-phase component; ( — ■) out-of-phase 
component of p^(r). 










Figure 27: Local work per cycle distributions and global works per cycle for the Rotor 67 
fan undergoing pure torsional vibrations about midchord at lo = 0.54 and a = 0 and 180 
deg (Np = 0 and 11). 
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Figure 28: Local work per cycle distributions and global works per cycle for the Rotor 67 
fan undergoing pure bending vibrations at u — 0.54 anc a — 0 and 180 deg (No = 0 and 
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